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The instability of uniform viscous flow under 
rollers and spreaders 


By J. R. A. PEARSON 


Imperial Chemical Industries Limited, Akers Research Laboratories, 
The Frythe, Welwyn, Herts.+ 


(Received 24 July 1959) 


When a thin film of viscous fluid is produced by passing it through a small gap 
between a roller or spreader and a flat plate, it often presents a waved, or ribbed, 
surface. An analysis is given here in terms of lubrication theory to show why in 
many cases flow leading to a uniform film is unstable. Account is taken of surface 
tension which proves to be a stabilizing factor. The most unstable values of the 
} wave-number, » (characterizing the disturbance), are calculated as functions of 
the dimensionless variable 7'/~U), and of the geometry of the system; 7’ is the 
surface tension, / the viscosity and Uy a representative velocity of the fluid. For 
the particular case of a spreader in the form of a wide-angled wedge, these 
predictions are compared with experimental observations. Agreement is obtained 
for values of 7'/uU, between about 10 and 0-1, but for smaller values of 7'/u~U 
it is clear that other considerations, involving only viscous and pressure forces, 
determine the nature of the secondary flow. 


1. Introduction 

Attempts to produce uniform thin layers of very viscous fluid by rolling or 
spreading often lead to an uneven or ribbed surface. Reference to this pheno- 
menon has been made in a variety of contexts; reports have come from the tin- 
plate (Chalmers & Hoare 1941), the printing ink (Sjodahl 1951), the photographie, 
leathercloth and paint industries (private communications), while the matter has 
also been discussed in the correspondence columns of the New Scientist (Letters, 
Nov., Dec., 1958). 

The average reader will perhaps be most familar with the case of ‘brush marks’, 
H which occur when thick paint is applied to a flat surface by means of a bristle 
brush. The resultant film, especially when just formed, usually presents a waved 
or ribbed surface where the lines of crests and troughs run roughly parallel to the 
direction of motion of the brush (see figure 1a, plate 1). A similar ribbing pheno- 
menon can be observed when films of highly viscous fluid are formed by passing 
the fluid through a small gap under large solid rollers; if the roller is moved 
perpendicular to its axis over a flat plate, without rotation, the operation is 
termed spreading; if the roller is rolled over the surface, the operation is termed 
rolling (see figure 1b,c). A diagrammatic indication of the method used for 
spreading and of the nature of the resulting film is given in figure 2. A more exact 


+ Now at The Metal Box Co. Ltd., Technical Engineering Division, Borehamwood, Herts. 
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indication of the shape of the meniscus in the region where the free surface forms 
is given also diagrammatically in figure 3 (see also figure 1e). 

Early quantitative experiments were carried out by A. J. G. Shaw (private 
communication) of the Research Department, Paints Division, Imperial Chemical 
Industries Limited, on a series on viscous liquids. Using circular cylindical rollers. 


Circular roller 
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FicuRE 2. The ribbing of a spread film. 
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Figure 3. Formation of free surface when ribbing occurs. 


of radii 3, } and 1}in., he produced layers of thicknesses varying from 0-0005 to 
0-016in. both by rolling and spreading. The viscosities of the liquids varied from 
about 5 to 300 P and the velocities of spreading and rolling were of the order of 
1—-20in./sec. In all cases of spreading and for relatively slow rolling he observed 
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that the emergent film was covered by regularly spaced lines of crests and troughs 
running parallel to the direction of motion of the roller. He was able to measure 
a characteristic wavelength for each set of conditions. In the case of rapid rolling 
the pattern became more confused with lines lying at varying small angles across 
the direction of motion. He found that for any given spreader (or roller) and given 
minimum gap width—in other words for any given geometrical configuration— 
these patterns were remarkably reproducible. The line spacing (wavelength) 
remained roughly constant whatever the translational velocity, U,, of the spreader 
(or roller) provided that it was sufficiently large to produce well-defined ribbing: 
as a result, he did not measure this velocity specifically. The patterns seemed to 
depend very little on the fluid used, provided that it was sufficiently viscous; there 
was, nevertheless, a small systematic variation with viscosity. Altering the 
minimum gap width or the radius of the roller did, however, have a marked effect 
on the wavelength and this dependence was measured for the range of values 
mentioned earlier. An example of the relation observed between line spacing— 
expressed as the number of crests per inch, n), and minimum gap thickness, hp, 
for a roller of ?in. radius is shown in figure 4. The curves for three liquids both 
spread and rolled are given. 

Other experiments were also carried out by E. Pitts and J. Greiller (private 
communication) of the Research Department, Kodak Limited on a pair of 
rotating rollers separated by a small gap and partially immersed in a bath of 
viscous fluid. In this arrangement the fluid was carried out of the bath between 
the rollers, and then split into a layer on each roller, thus being returned to the 
bath as the roller rotated. In this way a steady-state system could be achieved. 
They found that for sufficiently low rates of rotation a uniform two-dimensional 
flow was obtained, and that as the rate of rotation of the rollers was increased, 
a critical value was reached at which the flow became unstable and a ribbed 
structure similar to that described above was observed. The wavelength of the 
wave structure obtained was found to decrease—apparently in a discrete series 
of jumps—as the rate of rotation was further increased. The peripheral roller 
speeds used by Pitts and Greiller were in general lower than any used by Shaw. 

Dimensional considerations suggest that in the experiments described above. 
the various forces acting can be characterized in the following way: 


viscous forces = 


( 
inertia forces = O(pU?/hy). 
gravity forces = O(pq), 

( 


surface tension forces = O(T'/h@), 


where p is the density of the fluid, g is the gravitational acceleration, 7’ is the 
surface tension of the fluid and U,, hy and are as defined earlier. Ifh, is sufficiently 
small then viscous and surface tension forces dominate, and if further U, is large 
compared with 7 then viscous forces alone dominate. The characteristics common 
to all the industrial observations on ribbing mentioned earlier were that the gaps, 
ho, were very small and that the fluids were very viscous. The experiments of 


Shaw were such that viscous forces dominated everywhere, while those of Pitts 
31-2 
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and Greiller were such that viscous and surface tension forces were of comparable 
effect. 

This suggests that a theoretical treatment should be based on the equations of 
slow viscous motion, neglecting inertia and gravity forces, and that the effects of 
surface tension should be included in the boundary conditions. 


300 P 
80 P 
14P 


ho 


300 P 
14P 


(6) 


FIGURE 4. 7, the number of crests per inch, as a function of ho, the minimum gap thick- 
ness in 55 in., using a j in. radius roller. (a) Spreading. (b) Rolling. 


We present here an analysis based on the hydrodynamic theory of lubrication, 
which is an approximate theory of slow viscous motion applicable to the flows we 
are considering. The effect of surface tension on the free-surface boundary condi- 
tion has to be introduced in a simplified, semi-empirical form, because of the 
limitations inherent in the lubrication equations. 
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We look first for a possible two-dimensional steady-state solution, that is a 
solution corresponding to a uniform plane emergent sheet. We find that the 
equations and boundary conditions we have chosen lead to a degree of arbitrari- 
ness in this solution. Assuming, however, that one of this class of solutions will 
be relevant, we then study the stability of such flows to three-dimensional distur- 
bances of the type observed in practice. We find that the stability of the system 
is a function both of its geometry and of the single dimensionless parameter 
T |uU,; surface-tension forces exert a stabilizing influence on an otherwise wholly 
unstable flow. This general analysis is given in §2. We obtain expressions for the 
amplification factor, s, as a function of wave-number, n, of the postulated 
disturbance. 


he 


Emergent layer 


Flat glass plate 


Ficure 5. The wide-angled wedge spreader. 


In order to obtain solutions in terms of tabulated functions, a particularly 
simple geometrical arrangement is treated in detail in §3. This consists of a wide- 
angled wedge spreader, illustrated diagrammatically in figure 5. For any specific 
value of 7'/U,, and a suitably chosen value of b, (the point of furthest penetration 
of the meniscus), the value of n that maximizes s can be calculated. We assume 
that this is the value of n that would be observed in practice. 

)xperiments have been performed in these laboratories using a wedge-shaped 
spreader to compare the results obtained experimentally for with those predicted 
on the basis of the theory outlined above. These also are described in §3. The 
agreement is good for values of T'/U, between 0-1 and 10. 

It is concluded in §4 that the elementary analysis we have used involving 
surface tension is adequate to describe the ribbing that occurs when 7'/zU, is 
appreciable, but that other viscous effects, not included in the lubrication theory, 
would have to be considered in order to describe the secondary flow that takes 
place in the limiting case of 7'/U, > 0. 


2. General analysis 

The conditions encountered in flow under a roller or spreader, where the fluid 
is supposed incompressible, are precisely those to which the theory of lubrication 
applies (see, for example, Lamb 1932, Ch. XI). In this approximation, designed 
to deal with the motion of viscous fluid contained in the narrow regions between 
nearly parallel moving surfaces, inertia and acceleration terms are neglected, 
while pressure variations in directions at right angles to the moving surfaces are 


+ Pitts (1959) has developed an approximate analysis for predicting the point to which 
the meniscus will penetrate as a function of the dimensionless group T/wU,. 
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assumed to be unimportant. This means that a two-dimensional flow (of the type 
shown in figure 6) can be represented by a one-dimensional model, while the more 
complicated three-dimensional flow that results from the ribbing considered in § | 
can be represented by a two-dimensional model. This approximation can be 
regarded as adequate except in the regions where a meniscus forms, e.g. the line 
x = bin figure 6. This is because the parabolic type of flow which is relevant for 
regions far from the meniscus, where the gap is completely filled with fluid. 
changes sharply, in the region close to the meniscus, into uniform flow (with a free 
surface) along the two separating surfaces. The latter flow also satisfies the 
lubrication equations, but in a trivial sense only. The length (in the x-direction) 
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Moving surfaces 


\ \ 


Fluid | 
a 6 


FIGURE 6. Section of flow in the narrow gap between two arbitrarily moving surfaces. 


of the region directly affected by the meniscus will be of the order of the gap width 
at that point and is therefore small compared with other lengths involved, such as 
4 or the line spacing on the emergent sheet.t 

In this analysis we shall attempt to represent the whole effect of the meniscus 
formed at x= 06 by boundary conditions to be imposed on the lubrication 
equations at what we shall term a free boundary; this free boundary will for 
simplicity also be taken to be at 2 = 6. Our mathematical model (see figure 7) will 
therefore consist of two regions (neglecting for the moment the presence of a 
meniscus at 2 = —a): (i) x < b, where the gap is full of liquid and where the 
lubrication equations apply, and (ii) x > 6, consisting of two fluid layers of 
constant width moving with uniform velocity separated by an air gap. 

Although changing discontinuously the flows in these two regions must be 
matched. This is the crux of the problem: we must choose two suitable boundary 
conditions for the lubrication equations. There must be a mass-flux boundary 
condition which equates the amount of fluid reaching the free boundary from 
region (i) to the amount of fluid leaving the free boundary into region (ii). There 
must also be a pressure boundary condition, the pressure in region (i) at the free 
boundary being balanced by the sum of the (constant) pressure in region (ii) and 
the effect of surface tension forces acting on the free boundary. The particular con- 
ditions adopted in the following analysis will be discussed when they are chosen. 


+ This remark is based on empirical observations. 
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2.1. Steady two-dimensional flow 
To preserve the utmost generality in this treatment we consider the arrangement 
shown diagrammatically in figure 6. The two moving cylindrical surfaces are 
given by y = h,(x) and y = h(x), where hy(x) and h3(x) are everywhere small com- 
pared with unity. (The primes denote differentiation with respect to x.) The 
tangential velocities of the two surfaces are given by U, and MU, as indicated 


| | 
| 
Free boundary 


FIGURE 7. Mathematical model of flow shown in figure 6. 


(0 < M < 1). Using the usual assumptions of lubrication theory, we find that in 
the steady state the velocity in the x-direction must be of the form 
y) = U —E(x)}, (1 
where E(a) remains to be determined. The total mass flux, U(a), between the two 
surfaces for unit width in the z-direction (perpendicular to the z- and y-axes) at 
a section x is given by 
h(x) 
Ue) = = + M + (2) 
By continuity this must be constant, and if t, is the ultimate film thickness as 
x — 00 on surface | and ¢, the thickness on surface 2, then 
= U[t, + Mt]. (3) 
This is our first boundary condition. (¢, and ft, are the film thicknesses used for 
region (ii) in our model.) 
The equation of motion is given by 
dp 
de ~~ (hale) — hye) 
where pis the pressure and yz the dynamic viscosity. If we assume that the pressure 
in region (i) at the free boundaries is known, then equation (4) may be integrated 
using (2) and (3) to give 
b 
p(b)—p(—a) | (1+) 2(t, + | (5) 


—hy(x)}? — hy (x)} 
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Unless further information is provided a large degree of indeterminacy remains, 
for a, b, t,, tz, are all as yet arbitrary and are related by only one equation. 

One reasonable simplification is to suppose that a > oo and that p(—a) > 0. 
The pressure outside the liquid is taken to be zero. This corresponds to the case 
where the surfaces are emerging from a large reservoir of fluid, and is in fact 
reasonably representative of all the cases that have been examined above. 

But the indeterminacy still remaining in the specification of the section at which 
the free boundary—or meniscus—forms cannot be so simply resolved. The 
obvious refinement to attempt would be to replace the lubrication equations (1)- 
(4) by the Stokes slow-motion equations, and seek solutions of these using exact 
free surface boundary conditions on both velocity and stress. Using a stream 
function, this leads to the biharmonic equation. Unfortunately, the equation for 
the surface remains an unknown variable, and it is by no means certain that the 
solution of this more complicated equation would remove the indeterminacy. 

Pitts (1959) has described an approximate method for resolving this in- 
determinacy that he has applied with some success to the symmetric case of two 
circular rollers, where h,(a) = —h,(v) and M = 1. In this he makes use of the 
experimentally observed fact that the meniscus is parabolic near its intersection 
with the axis of symmetry; this equivalent parabola he characterizes by r, the 
length of its latus rectum. He then looks for that value of b (see figure 6) which for 
suitable r makes the lubrication solution (1) (written in terms of a stream function) 
satisfy the exact stress boundary conditions on the parabolic surface in the region 
close to the axis of symmetry. The predicted values of b (for values of 7'/uU, 
sufficiently high for the two-dimensional flow to be stable) were compared with 
experimental observations and fairly encouraging results obtained. Unfortu- 
nately, it is not clear how such a method could be applied to the asymmetric case 
of spreading, and this will not be attempted in $3. 

Hopkins (1957) has suggested that the meniscus will form at the first stagnation 
point (in the case of a sheet moving between two rollers); this seems to be con- 
firmed for the case 7'/jl), — 0 both by experimental results and by the approxi- 
mate theory of Pitts. 

Experimentally, the value of } is found to decrease as U, is increased (for fixed 
T and ) and tends to an asymptotic value, b,,. This general result is to be expected 
even on the basis of our approximate treatment, because the lubrication solution 
yields a region of reverse flow for all b greater than some value b,, b, being defined 
by the geometry alone; we would therefore expect the meniscus to be drawn 
towards the point },, the more so as /U) increases with respect to 7’, and indeed 
it has been roughly verified that b, = b,. 

One further remark can be made about the variables in (5) for the particular 
arrangements that have been studied. For rolling, WV = 1 and by symmetry we 
could put t, = f,, while for spreading, M = 0 and f, = 0. Thus if we fiz b, and 
assume that p(b) is known, then ¢, is given by (5) in terms of known quantities. 
and E(2) in (2) can be calculated. In this way, specification of 6 (assuming 
a -> 2) leads to a complete solution for the velocity u(x, y) in terms of lubrication 
theory. Little has been said about p(b) so far: the simplest assumption is to put it 
equal to zero. However, allowance can be made for surface tension effects as will 
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be seen in later sections, and indeed the inclusion of surface tension effects will 
be critical to the stability analysis that we shall carry out in the following 
subsection. 


2.2. The growth of small lateral disturbances 
We now superpose on the given basic two-dimensional steady flow infinitesimal 
lateral disturbances that are chosen to correspond to those described in §1. 
Since the equations describing the motion are linear there will be no interaction 
between the two motions in the main body of fluid. Hence it is the particular choice 
of boundary conditions that will be of critical importance. 
We therefore now consider a non-steady velocity field given by 


(6) 


[ha(x) — (x) 


where w and w are velocities in the x- and z-directions, respectively. There will be 
associated mass fluxes 


U(a2,2,t) = | udy = 4U,[ho(x) —h,(x)][1+ M +3{K(z, z, t)}] (8) 
h,(x) 
and W(z,z,t) = wdy = LU —h,(x)] J (x. z, t). (9) 
J hi(x) 
Continuity imposes the condition 
(10) 
Cx 02 
whence = 0. (11) 
while the equations of motion are given by 
Op dp dw 


We now characterize the disturbance by supposing the free boundary to form 


at the points b = b,+ce“cos nz, (13) 


where c < 1 and » is arbitrary. We suppose in consequence that the mass flows 
be given by K(x,z,t) = E(x) —cF (x) et cos nz, (14) 
J (x, 2, t) = cG(x) sin nz, (15) 
where E(x) is defined by (2), (3) and (5). By substituting (14) and (15) into (11) 
and using (2) and (3), we obtain from the continuity relation 


h,) F} = nG(h,—hy); (16) 
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by using the equations (12) and eliminating p we get 


(G 
where h=h,—h,. Finally. we obtain the following second-order linear dif- 
ferential equation for G(z). 
h’ 2h" 
—— G' —(n?+——]G = 0. 18 
+5 ) G=0 (18) 
We must now choose two boundary conditions to be imposed at the free boundary 
given by (13). 
The first boundary condition will be obtained by considering the force balance 
across the free boundary. We can without loss of generality put p = 0in region (ii). 
If we consider the two-dimensional case represented in figure 7, then the relevant 


relation is given by p(b)h(b) = —2T7. (19) 


where p(b) refers now to the pressure in region (i). This relation is consistent with 
the lubrication hypothesis in that inertia terms are neglected, and with the con- 
cept of a free boundary, thus avoiding any viscous components of stress.+ The 
three-dimensional case must include surface tension components that act in the 
xz-plane and by using the simplest possible representation we obtain 


p(b)h(b) = — (20) 


If we now linearize in c this pressure boundary condition (20) becomes 


2 & & | 
p(b) Tia +n (1 h 4 e cos (21) 
We may further assume that 
p(—a)=0, G(-a)=0, F(-a)=0. (22) 


This is in keeping with the observation that the lateral variations in layer thick- 
ness of the emerging film have no effect on the flow of the fluid in front of the roller 
or spreader, and the further observation that the nature of the emerging film is 
independent of the conditions existing in front of the roller or spreader provided 
only that a is large enough (i.e. a > b). If we therefore integrate the first of the 
equations of motion (12), linearize in c, use relations (6), (14) and (21) and subtract 
out the steady component satisfying (5), we derive our first boundary condition 
on the time-dependent component. This is 


G(b,) = —n | sath (bo) + 2) (23) 


It is here that dependence on the parameter 7'/U, is introduced. 


+ This simplication was suggested originally by Sir Geoffrey Taylor (see Saffman & 
Taylor 1958). 
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A second boundary condition is obtained by relating the rate at which fluid 
reaches the free surface to the rate at which it is carried away, and to the rate at 
which the free surface advances. We write 

(b) 


m(b) = m,(b) = 


t¥(b) 
h(b)° 


where ty and ¢3 are the instantaneous thicknesses of the layers formed on surfaces | 
and 2 along lines z = const. (These must not, for the present, be confused with 
the values ¢, and t, considered in the previous section, since a basic unperturbed 
flow may be possible for one value of 6 only.) For any value of z, the amount of 
fluid instantaneously arriving at b is given by 


h(b) [1+ M + 3{E(b)} + Ac{F(b)} e* cos nz]: 
the amount instantaneously carried away by the surfaces | and 2 is given by 
Up[m,(b) + Mm,(b)] h(d); 


the rate at which the free surface advances is given by 
ob 
— = cos nz. 
ct 
Mass conservation therefore requires that 
ob 
(l1—m,—m,)h + M +44 + cos nz]—U)(m, + Mm,)h.t 


Linearizing and using the steady-state result given by (2) and (3), we get 
Uy 


as the second boundary condition. Here the function Z, and hence E£’, is taken 
to be defined uniquely by (2), (3) and (13) whatever m, and m,. Unfortunately, 
our whole analysis so far does not enable us to say anything definite about 
and m;(b,) though we can suppose m,(b,) and m,(b,) to be known. However, 
if we use (24) only to select the most unstable wave-number, i.e. that value of 
n that maximizes s, then the value of the constant term that remains unknown 
is not relevant.{ It is only when we use (24) to provide criteria for absolute 
stability that we need to evaluate this constant. 

If we take m, and m, to be given as functions of b. by equating ff and ff to t, 
and fy, for all b. then relation (24) becomes simply 


6{1 — my (bo) — Ma(bo)} 
This is equivalent to supposing that all basic unperturbed flows are initially 
possible. If, on the other hand, we take m}(b,) = m3(b)) = 0, which corresponds 


(24, i) 


+ This relation assumes that velocities in the z-direction are not transmitted across the 
free boundary. This approximation will be valid provided the wavelength 27/n is large 
compared with h. This has been the case in all our observations. 

t This supposes that m; and mz; are to first order in c independent of n. 
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to the supposition that the shape of the free surface in the xy-plane for all 6 will 
be everywhere geometrically similar, then (24) becomes 

Us Ald 1 > 


In the application of the theory made in §3 we shall use both of these criteria. 


3. The wedge-shaped spreader 

In order to study the consequences of the analysis given in the previous section, 
a particular case, that of a wide-angled wedge, has been worked out in detail and 
the results compared with those obtained experimentally. The relations (24, i) 
and (24,ii) have been calculated for various values of 7'/~U, and 6b, regarding 
s as a function of n; 4”, U, and 7' can be taken as adjustable parameters; b, is not 
prescribed by our elementary theory, but it has been measured experimentally 
for various values of 7'/wU). It must be remembered that the analysis of §2.2 isa 
perturbation analysis and strictly speaking refers only to the growth of small 
disturbances on a steady two-dimensional flow. If the initial disturbances are 
very small and are taken to include components at all wave-numbers, then it is 
reasonable to assume that those wave-numbers which maximize s for given b, 
and 7'/~U, will dominate the early stages of the instability. But it is not necessarily 
true that the ultimate steady secondary flow that is attained will be characterized 
by the same wave-number, although other examples of instability in hydro- 
dynamics suggest that it may be. Nevertheless, since we are unable to carry out 
a full analysis of the secondary flow, we are bound to carry the comparison as far 
as we can—on the supposition that agreement between perturbation theory and 
observation in a range strictly outside that to which perturbation theory applies 
is better than no agreement at all. Ideally we should restrict comparison to those 
values of 7'/~U, for which the flow is just unstable; unfortunately, this range, for 
the linear dimensions chosen in our experiments, corresponds to that region where 
gravity forces become appreciable as an additional stabilizing factor. In view of 
all these uncertainties, it is perhaps rather surprising that the agreement obtained 
was so good. 

3.1. Analysis 

A diagram of the wedge used is given in figure 5. The equation for the gap width, 
h(x), is given by h(x) = ho 25) 
The free boundary is supposed to form at x = 6. (This corresponds to the case 
h,(x) = 0, M = 0, t, = 0, in §2.1.) We write ¢, = t, A = t/hg, k = a/h. 

From (1) we have 


y 
ule, y) = PY (26) 
and from (2) and (3) 3194 — hla 
For steady-state flow, we suppose that a > «, p(a) > 0, and that p(d) is given by 
2T 


(28) 
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This is the assumption that we introduced earlier in §2.2. The pressure condition 


(5) becomes _ Te  _f A__ 1 9 2a] 


and thus £(b) can be written 
— 6(1 + kb)? + 12(1+kb)—-3+7 


2T a 
where t= rd (31) 


and will be used henceforth as the relevant dimensionless number characterizing 
the relative effect of viscous and surface tension forces. 
The basic equation for small perturbations, equation (18), becomes 


kx 
G’ = 0. 32 
(1+ ais 
For x > 0, we write 1+ kx = ky to give 
’G 1d@ 
i nG=0 > &), 


and for x < 0, we write 1—kx = k€ to give 
These equations have as general solutions 
G = 
= + (ng), 


(€>k). 


(33) 


where J, and K, are modified Bessel functions. If we apply boundary condition 
22), then as x > —co and €> 0, G must > 0. Hence A, = 0. Also at x = 0, 
G and G’ must be continuous. Hence 


This yields a unique relation between A, and B,. Remembering that we are 
supposing the free surface to form at the point b = b,+ce“cosnz, the pressure 
condition (23) yields the relation 


Po) +B,K, (7 + nb, )| 


—n] (1 —m(by nto) |]. (35) 


where E(by) is defined by (30) and m(b,) = }[£(b,)+3]. Thus A, and B, are 
completely specified in terms of by, 7 and n, using (34) and (35), giving 
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n| E(by) + + — ([n/k] + 


A,=- 


The equations (24, i) and (24, ii) become 


+nbo) +2K,( 7 +nbo)] 


6n[1 — m(b,)] 


6n{ 1 —m(bp) (37, ii) 


(36) 


where A, and B, are given by (36). 

Because of the arbitrariness in by, we are not able to present complete stability 
diagrams, that is curves for s as a function n and 7, from theoretical considerations 
alone. Thus we cannot predict in advance that the motion will always be stable 
forall7 > some 7,,,, and unstable for 7 < 7,,,,,,even though experimental observa- 
tions suggest that this is true. Our mathematical model must use empirical 
information if it is to be compared with experimental results. 

Some typical examples of the curves obtained from (37, i) and (37, ii) are given 
by figure 8. Three separate values of 7 have been chosen; for each of these, values 
of b, which correspond to those observed experimentally (for particular values of 
a and h,) have been used. Figure 9 attempts to show how n(smax) (where (max) 
refers to that value of n that maximizes s) varies with b, for fixed 7: it is a fairly 
typical curve and shows how little n(smax) varies over quite a wide range of b, 
which includes the observed value. 

It can be seen that for all values of by, the hypothesis (i) leads to a more unstable 
situation than hypothesis (ii); this conclusion is fairly obvious on physical grounds 
since the second hypothesis supposes that increasing b, increases ¢ much more 
rapidly than does the first. 

The behaviour of s in the limit n = 0 is of some interest: the relevant relations 
that follow from (37,1) and (37. ii) are 


(bo) +7] 
6[1 — m(b)][(1/k) + Bo] 
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(38,i), and hence our mathematical model, implies that for all 7, there exists a 
b.,i, Such that s > 0 for all by < b.,,, while (38, ii) implies a range of instability 
that is contained within that obtained using (38,i). It can be argued that our 
supposed steady-state two-dimensional flow should never be characterized by a b, 
lying within the range of instability for zero wave-number disturbances, and 
indeed it was found that the observed values of by always satisfied this condition. 


b, observed 


FIGURE 9. 7(8max)/k as a function of by for 7 = 0-1. 


3.2. Experimental results 


For the purposes of experiment a compromise had to be made between values 
of h, and @ sufficiently large to allow of accurate machining and sufficiently small 
to ensure the validity of the lubrication equations and of our free boundary 
hypothesis. With this in mind, the values hy = 0-004 in. and « = 1/20 were chosen. 
Furthermore, the spreader could only be of limited size. The wedge was termi- 
nated at x = lin. anda = — 1 in. (using the notation of figure 5) by vertical sides, 
and was some 6in. wide in the z-direction. It was drawn, at measured constant 
speed, across a plate-glass sheet by a geared-down electric motor, a side flange 
ensuring that its motion was perpendicular to planes x = const. Instantaneous 
flash photographs were taken from underneath the plate to measure bo, an 
example of which is given in figure le; photographs taken from above by trans- 
mitted light were used to measure the line spacing, and an example of these is 
given in figure 1d. 

)xperience with circular rollers of different lengths suggested that the limita- 
tion on size in the z-direction would not have any appreciable effect, while the 
truncation of the wedge at x = — 1 was only expected to have noticeable effect at 
the lower speeds. This can be seen by solving equation (32) using a non-infinite 

value for a in (22); the effect is not large provided n is not small. The patterns of 
lines obtained, although not absolutely regular on any photograph, were repro- 
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Figure (plate 1). (a) Brush marks—bristle brush. (b) Spread film — circular roller. 
(¢) Rolled film— circular roller. (d) Spread film—wedge spreader. (e) Under view of wedge 
spreader during experiment. 
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ducible and yielded satisfactory estimates for n in each case. Three fluids were 
used; they are tabulated with their relevant physical properties below. 


Surface 
Viscosity tension 
Liquid (poise at 23°C) (e.g.s. units) 
L.C.I. silicone oil F110/1000 14:5 21 
Analar glycerol 11-6 63 
B.D.H. glycerol diluted with water 2-7 57 
50 r- 
/ 
/ 
H 
N 
20 
: 
1 + 
: 
10 
10 03 O1 0:03 0-01 0-003 


FicurE 10. N, the number of crests per inch, as a function of 7—comparison between 
calculated range of instability and observed values. —-——, calculated neutral stability 
curve; , calculated n(smax); L, glycerol; 3P; 1, glycerol 12P; T, silicone oil 15P. 


The combined results are presented in figure 10 were 1, )cerveq) 1S plotted as a 
function of 7. The values of Do (o4serveq) ae plotted as a function of 7 in figure 11. 

Owing to the large variation of viscosity with temperature and with water 
content of glycerol, and to the difficulty of preventing surface contamination and 
hence a reduction in effective surface tension, the calculated values for 7 are 
subject to a rather large margin of error, possibly even 20° in some cases. This 
does not mean that quantities were not measured accurately or that reasonable 
precautions were not taken where necessary, but it allows for the occasions when 
conditions happened to change between measurements of relevant physical 
quantities. Indeed, it is worth pointing out that extreme care was taken to polish 
the glass plate and to clean the metal surfaces on the wedge before each experiment 
and that fresh glycerol was used every time. 


3.3. Comparison of theory and experiment 


Because the analysis of §2.1 does not prescribe a value of ), for any particular 
choice of 7, we have used the observed values of b, as a starting point for the 
32 Fluid Mech. 7 
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perturbation calculations. Figure 11 shows the results obtained and compares 
them with the values that would have been predicted supposing the free surface 
to form at the point where reverse flow is just about to take place, namely where 
E(b,) = —1. It is seen that as 7 decreases the two curves tend to coincide. This is 
in agreement with the observations of Pitts and Greiller on two revolving rollers 
and is the analogue of Hopkins’s postulate. It can equally be verified that b, > b,, 
(see §3.1) in all cases. 


0-4 
03 
= 
E(6))=-1 
l 
O1 0-01 
T 


Figure 11. by) vs r-—observed values. 


The results of the calculations based on (38, i) and (38, ii) are shown in figure 8, 
for three particular values of 7 (7 = 1, 0-1 and 0-01) and appropriate values of bp. 
These give an indication of the large range of wave-number that is likely to be 
unstable. A more comprehensive indication of the predictions is given in figure 10 
where ”(Smax) and the entire range of unstable n are plotted as functions of 7. The 
width of the band is a function of the form assumed for m’(b) in (37) and in this 
figure the means of values obtained by hypotheses (i) and (ii) have been plotted. 
Also shown on this figure are the observed values of n, mentioned in the previous 
section. The agreement over a wide range of 7 is encouraging. It will be noticed 
that neither the calculations nor the experimental observations are given beyond 
the point 7 = 1. From the experimental point of view it was clear that a critical 
value did exist for 7(7 ~ 3) above which the two-dimensional flow was stable; 
however, the amplitude of the ribbing for values of 7 just below this critical value 
was so small that measurement of wave-number became extremely difficult; this 
also made it difficult to determine the critical point exactly. Also, in our pre- 
liminary discussion of various forces that might be relevant, we showed by 
dimensional considerations that gravity forces would be small compared with 
viscous and/or surface tension forces, provided hy were small enough. This result, 
like many dimensional arguments, is only partially true. If U, becomes sufficiently 
small, then the relevant length scale is not hy but h(by), and the ratio between 
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viscous and gravity forces is given by “U,/h?og and between surface tension and 
gravity forces by 7'/h?pg. In our particular observations, when 

T=1, T/h?pg = 10, 
and gravity forces are thus of the same order of magnitude as viscous forces. To 
emphasize this fact we have deliberately truncated figure 10 at 7 = 1. 

As 7 tends to zero, the observed value of n is seen to tend to a constant value. 
Because of the extremely high value for viscosity, or the extremely high value for 
velocity involved at low values of 7, reproducible steady flows were difficult to 
obtain with the given length of glass plate. In fact, b) was observed to vary 
between experiments; interestingly enough the values recorded for n varied 
simultaneously in the manner suggested by figure 11. 


4. Discussion 

The comparison made in §3.3 for the special case of the wedge-shaped 
spreader suggests strongly that a reasonable theoretical explanation has been 
given for the formation of the observed line structure where 7 > 0-03. The flow 
can be interpreted primarily as a balance between pressure and viscous forces, 
with surface tension exerting, through the boundary conditions, a stabilizing 
influence. In principle, exactly the same analysis can be carried out for arbitrary 
functions h,(x) and h,(x) and for arbitrary values of 1. However, even in the case 
of circular rollers, the complete solution of equation (18) does not seem readily 
expressible in terms of tabulated functions, and a complete numerical solution 
would be laborious. It is, nevertheless, reasonable to suppose that the results of 
Pitts and Greiller, where 7'/~U, = O(1), could be explained using such an analysis. 
The results of Shaw are to be interpreted as corresponding to the region 7 > 0. 
It is clear that for this range of values of 7'/~U,, the analysis is unable to represent 
the secondary flow that occurs; the predicted values of n(8max) tend to infinity and 
no longer correspond to observed values. This aspect of the problem was discussed 
in an earlier investigation (Pearson 1957), where it was concluded that the flow 
pattern in the columns of fluid (figure 3) exerts a dominant influence on the line 
spacing. 

From the purely practical (industrial) point of view, the relevant question to 
be answered is whether conditions can be so chosen that a stable two-dimensional 
flow is obtained even for small 7'/uU,. (The trivial answer that a stable flow can be 
obtained for sufficiently small U, is of little assistance.) In most contexts the 
problem will be tackled and probably solved (or, if not, circumvented) by 
methods of trial and error, since the possible range of variation of physical 
circumstances will be strictly limited in any one case. 

One method that proves efficient was mentioned in an earlier investigation 
(Pearson 1957). The roller or wedge is replaced by a flat plate ending in a knife 
edge. This plate is inclined at a small angle to the horizontal such that h’(x) is 
negative everywhere, and such that fluid emerges at the point of minimum gap 
width. Clearly the analysis given above cannot be applied to this case since by will 
coincide with the end of the plate and h’(z) will be discontinuous at by. On the 
other hand, a simple physical argument in terms of pressure can be given to 
explain why the flow should be unstable in some cases and stable in others: 
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consider the flow shown in figure 7 and let 6, represent the position of the free 
boundary in a steady two-dimensional flow; if the pressure gradient dp/dx within 
the fluid is negative near the point x = b, in the steady two-dimensional case 
then the result of part of the free boundary moving up to the plane 
x=b,+6b (db > 0) 

is that the induced pressure distribution forces fluid laterally towards these 
points. In other words, it tends to provide just the excess of fluid that is required 
to maintain the disturbed boundary. Now, in the case of the roller or the wedge, 
dp/dz is negative and instability is observed. But where h’(b,) is negative, as we 
propose using a knife edge, dp/dz is positive and the flow is inherently stable. 

Another technique that has been employed is to make M = —1. This again 
puts the physical circumstances outside the scope of the analysis in §2 because 
there is now no clear distinction between an upstream and a downstream side. 
Nevertheless, by putting t; =¢, one can obtain a steady solution in which 
E(x) = Oand by following through the perturbation analysis, it can be shown that 
the flow is stable provided that d?/dz*(1/h) is negative at x = b). No comparison 
with observation has yet been made in this case. 

An incidental result of this work was that it suggested a new explanation for 
brush marks. In the past, it had commonly been supposed that these irregu- 
larities were due to the nature of the brush itself, formed as it is from a large 
number of bristles of roughly equal length, the whole assembly being relatively 
pliable. It now seems as if this is not so and that brush marks are merely a further 
manifestation of a phenomenon common to most spreading devices. 


Since much of the preliminary work on this problem was carried out in the 
Research Department of Imperial Chemical Industries Limited, Paints Division 
it is fitting that I should first acknowledge the stimulus given and information 
provided by Mr C. I. Snow, Mr N. D. P. Smith and Mr A. J. G. Shaw of Paints 
Division. 

Next I wish to express my indebtedness to Sir Geoffrey Taylor who has been 
kind enough to discuss the problem on numerous occasions and to make available 
many pieces of information that he has obtained on this and related matters in 
recent years. 

Finally, I wish to thank those members of these laboratories who have helped 
me with the experimental arrangements—in particular, Miss E. M. Gibbens, who 
assisted with measurements and calculations, Mr K. D. Cooper who was re- 
sponsible for most of the photographic arrangements, and the workshop staff 
who made the spreaders and rollers used here. 
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Instability of a viscous liquid of variable density 
in a vertical Hele-Shaw cell 


By R. A. WOODING 


Emmanuel College, Cambridge 
(Received 8 July 1959) 


Approximate equations of motion, continuity and mass transport are given for 
a viscous liquid of variable density moving very slowly between vertical and 
impermeable parallel planes. These equations are used to calculate approximate 
stability criteria when the liquid is at rest under a vertical density gradient. The 
results are applicable to the problem of the stability of a viscous liquid of variable 
density to two-dimensional disturbances in a porous medium. 

An exact stability analysis for the liquid between parallel planes is also given, 
and expansions in powers of the disturbance wave-number are obtained for the 
critical Rayleigh number at neutral stability. The previous approximate results 
are found to correspond to the leading terms of the series expansions. For the 
most unstable type of disturbance, the velocity distribution closely resembles 
plane Poiseuille flow, which was the form assumed in the approximate equations. 

An asymptotic expansion is derived for the critical Rayleigh number at neutral 
stability in a long vertical channel, or duct, the cross-section of which is a thin 
rectangle. The typical neutral disturbance possesses a ‘boundary layer’ at each 
end of the cell cross-section, and this has a small stabilizing effect. 

The critical Rayleigh number for a long vertical channel of rectangular cross- 
section is found experimentally by comparing the density gradient of the liquid 
in the channel at neutral stability with the corresponding density gradient in a 
vertical capillary tube. There is better agreement with the exact theory than 
with the approximate theory, the experimental result being about 4°, higher 
than the value predicted by the ‘exact’ asymptotic expression, and about 10% 
higher than the value predicted by the simple approximate theory. 


1. Introduction 

It is well known (see Lamb 1932, $330) that the flow of a viscous liquid at very 
low Reynolds numbers between parallel planes may be described by approximate 
equations of motion which have a simple form. Consider, in particular, vertical 
parallel planes which are separated by a small distance 2h. Take the origin of 
a rectangular co-ordinate system (X, Y, Z) midway between the planes, with OY 
normal to them and with OZ directed vertically upwards. It will be assumed, 
first, that the component of velocity normal to the planes vanishes, so that the 
pressure P is a function of X and Z alone, and secondly, that the components 
of velocity in the directions of OX and OZ are slowly varying functions of 
those co-ordinates. Then the velocity distribution is nearly parabolic in the 
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Y-co-ordinate, and if mean values are taken by integrating with respect to Y, the 
mean-velocity vector is everywhere parallel to the hydraulic head. The approxi- 


mate equations of motion are 
h? oP h? (oP 
where u and w are the components of the mean velocity, P, p and yw denote 
respectively the mean pressure, density and viscosity of the liquid, and g is the 
acceleration due to gravity. 
When p and are constants, equation (1) and the equation of continuity 


ou ow 
2 
ax (2) 


together show that the mean motion is irrotational and that a mean-velocity 
potential exists. Hele-Shaw (1898) has used this principle experimentally to 
render visible the streamlines of two-dimensional ideal potential flow about 
cylindrical bodies, and has treated, in a similar manner, lines of magnetic 
induction about elliptic cylinders and cylindrical shells (Hele-Shaw & Hay 1900). 

Saffman & Taylor (1958) have observed that the equations (1) are identical 
with Darcy’s law for two-dimensional flow of a viscous liquid through a porous 
medium of permeability 442. Using this analogy, they have made theoretical and 
experimental investigations of the stability of the interface between two im- 
miscible viscous fluids in motion through a porous medium in a gravitational 
field. 

Since the observation of flow phenomena within porous materials is difficult, 
the use of a Hele-Shaw analogue in certain experimental problems could offer 
attractive possibilities. The work described in this paper arose from the investiga- 
tion of the instability of a liquid of variable density in equilibrium under gravity 
in a vertical tube filled with porous material (Wooding 1959). The Hele-Shaw 
analogue is a long vertical channel, with a horizontal cross-section in the form of 
a very narrow rectangle. 

Exact calculations of the stability of a viscous liquid of variable density 
between parallel planes have been made by Ostrach (1955) and by Yih (1959). 
However, these authors did not discuss the most unstable forms of infinitesimal 
disturbance that can arise. 

A more realistic approach has been adopted by Bentwich (1959), who examined 
the inhibiting effect of horizontal magnetic fields upon the gravitational instability 
of fluid in a verticai tube of square or rectangular cross-section. Bentwich did not 
extend his work to the case of either infinite parallel planes or a Hele-Shaw cell. 
As these examples are of some importance, the gravitational stability theory is 
discussed in this paper. 


2. Approximate stability theory 


When the space between two parallel planes is filled by a liquid of variable 
density, the nature of the density distribution is influenced by both the mole- 
cular diffusion and by mechanical dispersion arising from the slow motion of the 
fluid. The relative importance of these effects can be estimated by adapting the 
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work of Taylor (1953, 1954a) and of Aris (1956) on mass transport of solute in 
tubes to the case of unidirectional flow between parallel planes. For simplicity, 
suppose that the fluid motion is in the Z-direction, and that the concentration of 
solute is a slowly varying function of Y and Z only. If gravity effects do not 
modify the velocity profile appreciably, it is found that the dissolved material is 
dispersed relative to a plane Z = wt, moving with the mean velocity w, by a 
process resembling molecular diffusion but with a diffusivity 


2 wh? 
105 « 
where « is the molecular diffusivity. This expression shows that, if the mean 
velocity is so small that the relation 


is satisfied, the effect of mechanical dispersion of the solute is negligible compared 
with molecular diffusion. When this condition applies, the density p can be 
treated as a slowly varying function of X and Z, and the approximate equations 
of motion and continuity, (1) and (2), will apply. A further equation, the equation 
of mass transport, is then found to be 


= taza)? (3) 
where u and w are the X- and Z-components of the mean velocity, as before. The 
equations (1), (2) and (3) are identical with the corresponding equations of two- 
dimensional motion in a porous medium of permeability 4h? and unit porosity, 
saturated by a liquid of effective diffusivity x and of variable density p. 

At this stage, the equations are seen to possess a shortcoming. Because 
inertial terms are neglected in (1), the system of equations is of first order in the 
time, and initial values of velocity and density cannot be assigned independently. 
Clearly, the neglected terms must be significant in the initial stages of the motion. 
However, it is easily shown that the inertial terms become negligible within a 
period of time of order h?/v seconds, where v is the kinematic viscosity, after which 
time the equations (1) become valid. 

Now consider a viscous liquid at rest between vertical parallel planes. For 
equilibrium, it is necessary that the unperturbed values of the pressure P and the 
density p should be functions of Z only. It will be assumed that the density 
gradient dp/dZ is a constant, and that the total variation in density is sufficiently 
small for p, and quantities depending upon p such as viscosity and diffusivity, to 
be considered constant to first order. 

It is convenient to introduce dimensionless spatial variables (x, z) = (X/h, Z/h) 
using h (equal to half the distance between the parallel planes) as the typical 
length unit. 

If a small perturbation of arbitrary form is added to the equilibrium value of 
the density, the pressure field will be perturbed and a small velocity (u,w), a 
function of x, z and t, will appear. The linearized equations governing these small 
disturbances are determined from (1), (2) and (3). Suppose that one Fourier 
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component of the density disturbance is of form cos ax cos yz e“, where « and y 
are dimensionless wave-numbers in the x- and z-directions, respectively. Then the 
linearized perturbation equations form a compatible system provided that 


h*w a*X/3 
K (4) 
_ dp gh* 
where A= dZ UK (5) 


is the Rayleigh number. 

Equation (4) is, of course, identical in form with the corresponding compatibility 
equation for the stability of a viscous liquid in a porous medium (Lapwood 1948; 
Wooding 1959). The exact correspondence follows when h is replaced by a length 
unit 6, determined by boundary dimensions in the (X, Z)-plane. Dimensionless 
wave-numbers «,, y, are defined with respect to b, giving a,/b = a/h, and 
similarly for y,. In place of (5), one has 


dp g(h?/3)b? Ab 


An = RP? 


where A, is defined as the Rayleigh number for a porous medium of permeability 
1h?, saturated with a liquid, of viscosity ~, which diffuses through the porous 
material with an effective diffusivity x. 

At neutral equilibrium (o = 0), equation (4) gives 


3(ax? 


A ; 


(6) 
This formula shows that, when y is held constant, the minimum value of A arises 
fora = y, corresponding to square convection cells. However, if «is held constant, 
A has a minimum value when y = 0, the fluid motion being in the form of long 
vertical columns. The first case arises when the vertical Hele-Shaw cell is of 
finite height, but of infinite horizontal extent—a situation which is analogous 
to that of Rayleigh instability in a porous medium (Lapwood 1948). The second 
case arises when the vertical cell is of finite width but of infinite height, the 
situation being analogous to that ofa long vertical tube filled with porous material 
(Wooding 1959). In each case, the boundary conditions along the edges of the 
strip must be expressed in an approximate form involving mean values of the 
density p and the velocity (w,w). At a horizontal edge, the vanishing of the 
horizontal velocity component w is ignored, since the effect of the non-slip condi- 
tion becomes inappreciable at a short distance from the edge. Similarly, the 
vanishing of w is ignored at a vertical edge. 

Some comparisons will be made later between the above approximate results 
and the results of ‘exact ‘ calculations in §§3 and 4. 


3. Exact theory of neutral stability 

The approximate stability theory, based upon equations (1), (2) and (3), is 
satisfactory provided that the wavelength of the disturbance in each of the X- 
and Z-directions is large compared with the spacing between the parallel planes. 
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At shorter wavelengths, the analysis breaks down because the disturbance is no 
longer of a quasi-two-dimensional character. The effects of shear stresses in the 
X- and Z-directions, and of molecular diffusion in the Y-direction, become 
apparent at short wavelengths, and one would expect to find the physical system 
rather more stable than the approximate theory predicts. 

Yih (1959) has given a proof, along the same lines as Pellew & Southwell 
(1940) for Rayleigh instability, of the ‘principle of the exchange of stabilities’ 
for a viscous liquid confined by vertical walls of unlimited height and heated from 
below, when the disturbance is periodic in the vertical direction, and also for walls 
of finite height bounded at the ends by horizontal conducting planes. In the 
latter case, periodicity in the vertical direction is not a necessary assumption. 
Apart from the above restrictions, Yih’s proof is quite general, and is applicable 
to the work which follows. Then, since the behaviour of all neutral or amplified 
disturbances is aperiodic with respect to time, one may impose the condition 
o/ct = 0 for neutral stability. 

Taking Cartesian axes as before, with the origin O midway between the parallel 
planes, OZ directed vertically upwards and OY normal to the planes, one has the 
following differential equations governing a neutrally stable disturbance 


divu = 0, (7) 
= 0, (8) 


where p is the primary density distribution, # and p are perturbations in the 
density and pressure, and u = (u,v, w) is the velocity perturbation. In equa- 
tion (8), & = (0,0, —g). 

The method of solution used here will follow that of Hales (1937), who examined 
certain problems in the stability of liquid in vertical tubes. 

To eliminate p, one takes the curl of (8), and obtains equations for the two 
non-zero components of the vorticity 


10 
(10) 
ox 
Elimination of u and v between (10) and (7) gives the usual result 
= (11) 


where V? = V?—0?/0Z?. Then, if the spatial differentiations are referred to non- 
dimensional variables (a, y,z) = (X/h, Y/h, Z/h), (9) and (11) can be combined to 
ve the familiar equation Vow = (12) 


and a similar equation for ®. Here A is the Rayleigh number defined in (5). 
At the vertical insulating walls, the boundary conditions give 


u=v=w=0 and ed/ey=0 when y= +1. (13) 
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Separable solutions of the equations, of suitable form, can be found by as- 
suming that a disturbance can be Fourier analysed into components with respect 
to both x and z. For example, let 

u = U(y)sin ax sin yz, ° 
v = V(y) cos az sin yz, | 


(14) 
w = W(y) cos ax cos yz, 


d = O(y) cos ax cos yz. 


These represent components of a disturbance of cellular form in the (a, z)-plane, 

having dimensionless wave-numbers « and y, respectively, The relative phases 

of u,v, wand &, in the x- and z-directions, are determined by the relations (10), ete. 
From (12) and (14), W(y) obeys the equation 


(D? —a?—y?)3 W = A(D?— a?) W, (15) 
where D = d/dy. (When one puts a = 0, this equation corresponds to the equa- 
tion (40) used by Yih (1959). The disturbance is then of two-dimensional type, 


and depends upon y and z only.) Yih has shown that solutions of equation (15) 
are of form exp(+;y), where 


(t= 1,2, 3), (16) 
m, being a root of the indicial equation 
m3 — Am — Ay? = 0. (17) 


The form of equation (15) is such as to entail separate treatment of even and 
odd solutions for W(y). If it can be assumed that the disturbance of Hele-Shaw 
type introduced in §2 leads to the greatest instability, it appears that the even 
solution 


3 
W(y) = & A;cosh py (18) 
i=l 


is of greatest physical interest. For the corresponding odd solution, the hyper- 
bolic cosines are replaced by hyperbolic sines. 

In order to impose the boundary conditions (13), it is necessary to solve (9) for , 
and (10) for U and V, in terms of the solution (18) for W. Use is also made of the 
equation of continuity (7) and the symmetrical nature of the boundary conditions. 
It is found that 


2 3 
ere 

dZ mM; 
ow 4 cosh py 

yU = cosh ; (19) 
-yV=A sinh ay — A; ) 
i=1 


mM; 


where a = (a?+ 2), and C is a constant. The boundary condition U = 0 on 


y = +1 gives A; cosh 


cosha; 
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If the boundary conditions V = W = D@ = 0 ony = +1 are now imposed, the 
resultant characteristic equation is of the form 


|\h Z| 
M, M, M, | =0, (20) 
im 
where L,; = p,tanh tanha [a,) 
M, = p?-2?, (21) 
and N; = m? p; tanh p,. | 


When the odd solution for W(y) is considered, the same characteristic equation 
is obtained, but with cosh p, and sinh p; interchanged, and with tanha replaced 
by cotha. When a tends to zero (in this latter case) tle characteristic equation 
reduces to that given by Yih (1959, equation (46)). While still discussing three- 
dimensional disturbances, Yih made the assumption that the horizontal com- 
ponent of velocity parallel to the boundary planes was identically zero. However, 
it is clear from (19) that U does not necessarily vanish. 

The method of Yih can be used to show that, if « is assumed to have a fixed 
positive value, the critical Rayleigh number is stationary with respect to varia- 
tions in y when the wavelength of the disturbance in the z-direction tends to 
infinity. If the condition ¢A/oy = 0 is imposed, it is easily shown, from (16) and 
(17), that both ep,/éy and ém,/cy (t = 1, 2,3) contain y asa factor. Then, from (21), 
each of 0L,/dy, 0M,/ey and ON,/cy contains y as a factor. Following Yih, one 
differentiates the characteristic equation (20) with respect to y and puts 
cA/oy = 0. Each of the three determinants so obtained contains a column 
proportional to (aL, 

oy ey} 
whence the sum of the three determinants must vanish for y = 0. Thus the 
characteristic equation (20), and the condition dA/éy = 0, can be satisfied 
simultaneously for y = 0. Since positive values of y involve larger gradients of 
velocity and density than when y = 0, it appears that the stationary point y = 0 


gives a minimum value of A. 
In the special case y = 0, the roots of the cubic equation (17) reduce to 


(¢ = I, 2, 2), 


m, = At, m,=A-+*, m,=0, 


and the equations subsequent to (17) are found to simplify considerably. 
Equations (19) and (18), together with the boundary conditions, give 


coshp  coshq 


where p = p, = (a2+A?)}, q = py = («2—A#)? and the equation (20) reduces to 
ptanh p+qtanhg = 0. (23) 


For small a, A may be expanded in a series of ascending powers of a. The leading 
term, say Ap, of the series is given by solutions of the equation 


Ai (tanh Ad — tan = 0, (24) 
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for which the root of smallest magnitude is zero. Then A is of order a?, and the 
transcendental functions in (23) may be expanded in series of powers of their 
arguments, leading to the result 


A = 3a7(1+ + O(a). (25) 
A comparison with the result of §2—that A = 3a? when y = 0—shows that the 
simple approximate analysis gives the first term of the series (25). 
It is of interest to note that, when both « and y are small but non-zero, an 
expansion of the characteristic equation (20) gives 


3 2 ay2\2 
2 


+ terms of fourth order in « and y, 
for the smallest root of A, which corresponds to the equation (6) calculated from 
simple theory. If « has a constant value, this expression shows at once that the 
smallest root has a minimum value at y = 0. 

An expansion of w in (22) for the case Ay = 0 and y = 0 leads to 


w = {—34a(1—y*) + O(a3)} cos ax, (26) 


that is, the velocity distribution is almost parabolic with respect to y. This 
supports the assumptions of §2 that the fluid motion resembles Hele-Shaw flow. 

For odd solutions, in the case y = 0, the hyperbolic cosines in (22) are replaced 
by hyperbolic sines, and the characteristic equation (23) is replaced by 


Pp q 


where p and q have the same meanings as before. If a series expansion in powers 
of ~? is assumed to exist, the leading term (A,) of the expansion is found to be 
given by solutions of 


0, (27) 


tanh Aé + tan 

Here the smallest root cannot be equal to zero, and it is, in fact, approximately 
equal to (0-75287)! 31-29. 

The numerical value A ~ 31-29 is given by Ostrach (1955) and by Yih (1959) as 
the critical value of the Rayleigh number for neutral stability. However, even 
when « = 0, this mode of disturbance is not as unstable as that arising when one 
takes the zero root of equation (24) with « finite but small. A physical explanation 
for the difference is readily found. Ostrach and Yih consider disturbance modes 
of two-dimensional form, the motion being restricted to the (y, z)-plane. Factors 
contributing to the stability are shear stresses arising from velocity gradients in 
the y-direction and molecular diffusion in the y-direction. These factors apply for 
the most unstable disturbance of this type, when the fluid moves in two columns, 
one rising, say, in the range 0 < y < 1, and the other descending in the range 
0 >y > —1. However, when the velocity disturbance is of Hele-Shaw type, the 
density perturbation is almost constant with respect to y; the physical factors 
contributing to the stability of the liquid are shear stresses induced by velocity 
gradients in the y-direction, and molecular diffusion in the z-direction. When the 
wave-number « tends to zero, the latter effect becomes negligible. 


= 0. (28) 
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An alternative treatment of the case y = 0 makes use of the fact that equa- 
tion (12) simplifies if w is independent of z. The method was used originally by 
Sir Geoffrey Taylor for the problem of stability in a vertical tube. In the simplified 
form, (12) becomes 

, 
This expression states that (V{—A) w is an harmonic function of x2 and y which 
vanishes at the boundaries (from (9) and (11) and the boundary condition w = 0) 
and therefore vanishes everywhere. Then the fourth-order equation 


(Vi-A) w =0 (29) 
is equivalent to (12). This is the same differential equation as that governing the 
transverse displacement of a freely vibrating thin elastic plate. 


When the boundaries are two vertical and parallel insulating planes, the 
boundary conditions (13) reduce to 


oy 
and solutions of equation (29) give the results described in equations (22) to (26). 


It is perhaps worth noting that, if the boundaries are two vertical and parallel 
conducting planes, the even solutions of (29), with boundary conditions 


w = — (V?w) = 0, 


w = Viw = 0, 
are of form w = COsry COS ax, 


where r2+a2 = At and r = (n+4)7 (n = 0,1,2,...). The most unstable case 
arises when n = 0, and 
A = (42? + > 6-09 


when « + 0. In this example, A tends to a finite limit because the effect of mole- 
cular diffusion to (and from) the conducting walls remains significant when a 
becomes very small. Again, the velocity distribution is approximately parabolic 
in the y-direction. Ostrach (1955) briefly examined the two-dimensional modes 
of disturbance between vertical parallel conducting walls in connexion with the 
thermal instability of a viscous fluid, but the above more unstable mode was not 
discussed. 


4. Instability in a long vertical channel of rectangular section 


Some consideration will now be given to the calculation of the criterion for 
neutral stability in a vertical Hele-Shaw cell which is closed by boundaries in the 
(x, z)-plane. Two basic boundary configurations have been mentioned briefly 
in §2. 

Suppose that vertical parallel planes, of width 2b and spaced a distance 2h 
form the sides of a Hele-Shaw cell of great depth, the cell being closed by im- 
permeable vertical boundaries at the sides x = +L (L = b/h> 1). The dimen- 
sionless wave-number y in the z-direction will be assumed to be zero, but the 
wave-number « in the x-direction must be finite. As L is large, the problem of 
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determining the Rayleigh number A at neutral stability is similar to the problem 
of determining the characteristic frequency of a narrow rectangular plate, 
vibrating freely in the lowest transverse mode for which the average displacement 
is identically zero. From previous considerations, it would appear that a distur- 
bance motion antisymmetric about « = 0, but symmetric about y = 0 (taking the 
origin of dimensionless co-ordinates midway between the planes, as before) will 
furnish the lowest eigenvalue for which the mean displacement of the plate, or the 
net vertical flux of the fluid, will vanish. 

The appropriate differential equation for the vertical velocity w is given by (29), 
and the boundary conditions upon the four vertical insulating walls are 


(30) 


»= —(V?w)=0 on 
u iw) 


This problem can be treated by Taylor’s classical method of expansion in 
Fourier series (Taylor 1933). Consider the expression 


L 


wn 


-coshp, coshq,, , 


sinhr,, x sinh x | 
y}, (31 
where a,, = (n+ $)7, and A, and B, are constants. Each term in (31) satisfies the 
boundary conditions = 0 on x = +Landw=0ony= +1, and the 
differential equation (29) is satisfied provided that 


a? /L?+ Ab, = a;,/L? — A, ) 


n~ 


(32) 


It is necessary to determine non-trivial values of the constants A, and B,, in such 
a waythat (31) satisfies the conditions w = 0on x = +L and (é/ey) Viw = 0 on 
y = +1. This will be possible only if a certain characteristic equation, which 
relates A and L, can be satisfied. The values of A, and 6, are determined as 
follows. Since w is an even function of y which vanishes on the boundaries 
y = +1, it is clear that w may be represented on the boundaries x = +L by a 
Fourier series in cos a,,y, where m = 0, 1, 2,.... The vanishing ofwonz = +L then 
requires that each coefficient of the Fourier series so obtained should be zero. 
Similarly, (¢/cy) Vw may be represented on the boundaries y = + 1 by a Fourier 
series in sina,,2/L. The vanishing of (¢/ey) Vw on y = +1 then requires that 
each Fourier coefficient be zero. Proceeding in this way. one obtains two infinite 
systems of equations 
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The coefficients of A, and B, are 


(a?,/L? + a2)? — 2’ 

= Das 


mn 
2 r 


L 


mn 


(34) 


> 
n 8 n 


nn being the Kronecker delta. 
Equations (33) are homogeneous and, for non-zero solutions for A,, and B,, it is 
necessary that the characteristic determinant should vanish: 


Poo 0 Pou 
S=| 1 Py (35) 

| 


where P.,,, = L,,,/K i N}.. Itis easily shown that the sum of the moduli of the 
non-diagonal terms form a convergent double series. It follows that the deter- 
minant (35) is convergent (Whittaker & Watson 1950, p. 37). 

Now, it can be shown that, when LZ > 1, the characteristic determinant in (35) 
can be evaluated to give the first few terms of the asymptotic expansion for A. 
When A represents the lowest eigenvalue, the final result is 


(1 12s 108s? 


where a) = $7 and 
& 10045. (37) 

Each term in the s-series (37) arises from a term in the B,,-series in (31). Now, 
the form of the B,,-series is such as to give an exponentially small contribution 
except in a region close to each of the ends (2 = +L) of the rectangular cross- 
section. It follows that the fluid motion near 2 = +L is of boundary-layer type, 
the thickness of the layer being of order h. It is to be concluded that the terms 
involving s in (36), which increase the critical value of A, are contributions due to 
these layers. As expected, the result of putting s = 0 in (36) gives an expression 
of the form (25), which is applicable when no boundary layer is present. 

The analogy of the above theory with that of a thin elastic plate vibrating 
transversely shows that the same end-phenomenon will appear in the latter 
problem. 


5. Experimental work 


The experimental method used by Taylor (19546) to measure the density 
gradient of a solution in a vertical tube under conditions near neutral stability, is 
easily adapted to the case of a vertical channel of rectangular section. 
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An experimental apparatus employing two tubes connected to a reservoir A is 
shown in figure 1. The purpose of this arrangement is to compare the density 
gradient at neutral stability in the rectangular channel B with the corresponding 
density gradient in the capillary tube C, for which the value of the critical 
Rayleigh number is well established. Initially, the vertical channel and the 
capillary tube are filled from the bottom with a transparent liquid of uniform 
density p). A coloured solution of slightly greater density p, is poured into the 
reservoir A. Since the system is dynamically unstable, convection currents are 


B c 
A 
A 
FIGURE 1 FIGURE 2 


Ficure 1. Apparatus used for the comparison of density gradients at neutral stability. 


A = 100ml. reservoir, B = rectangular channel, C = capillary tube. 
Figure 2. Cross-section through the rectangular channel. A = glass strips, B = 30s.w.g. 
wire, C = cement. 


set up in both vertical tubes. When equilibrium has been reached, the density of 
the liquid in each tube decreases with depth at very nearly the critical rate, from 
the value p, at the reservoir to the value p, at an appropriate point some distance 
down each tube. Below those points, the liquid has uniform density py. Obviously, 
the depth of penetration of the colouring material indicates approximately the 
position at which the density of the liquid has decreased to py. If the dimensions 
of the channel and capillary tube are known, the ratio of the respective Rayleigh 
numbers can be calculated from the formula 


A hs 
A, 4a 


a 


(38) 


where A is the critical Rayleigh number in the vertical rectangular channel, 
2h is the spacing of the parallel planes, and Z,, is the depth of penetration of the 
upper liquid into the channel at neutral stability, A, ~ 67-94 is the critical 
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Rayleigh number in the vertical capillary tube (from theory, Taylor 19545), 
a is the radius of the tube and Z, is the depth of penetration of the upper liquid 
into the tube at neutral stability. An advantage of this comparison method is 
that the density, viscosity and diffusivity of the liquid need not be known when 
calculating the ratio A/A,. 

The rectangular channel used consisted of two flat glass strips about 40cm 
long, which were cemented together with Araldite (a commercial two-component 
cernent) using two stretched 30s.w.g. wires of diameter 0-0315¢em as spacers 
(figure 2). This was not an ideal arrangement, as the spacing between the wires 
was not quite constant at all points along the channel. However, the variation in 
spacing was not great, and it was sufficient to take an average value in the 
numerical calculations. A further disadvantage arose from the fact that the ends 
of the internal cross-section were not truly rectangular, but, as the spacing 
between the two wires was more than twenty times their diameter, the effect of 
the shape of the ends upon the stability characteristics should be small. 

The average width of the internal cross-section of the channel, shown in 
figure 2, was obtained by measurement of enlarged photographs; the spacing of 
the glass strips was obtained from micrometer measurements. In the case of the 
capillary tube, the internal volume was found by filling the tube with a measured 
volume of water, and from this the average radius was calculated. The results 
were 2Lh = 0-675 cm for the width and 2h = 0-0320cm for the thickness of the 
cell, and a = 0-165 cm for the radius of the capillary tube. 

Since the critical Rayleigh number for this cell in the vertical position was of 
order 10-2, a fairly viscous liquid was chosen—a water-glycerol mixture in equal 
proportions by weight, having a viscosity of about 0-06 poise. Portion of this 
mixture was coloured with methylene blue for addition to the reservoir at the top 
of the apparatus. The density ratio between the coloured and clear solutions was 
adjusted to about 1-002 by dissolving a small quantity of sodium chloride in a 
further portion of the clear mixture, and adding it drop-by-drop to each solution 
in turn. 

During the progress of the experiment, the apparatus was kept in a darkened 
room in which temperature variations were small. A period of about 12 days 
elapsed before convective motion ceased in the rectangular channel. 

After equilibrium had been reached in both the capillary tube and the rect- 
angular channel, the respective depths of penetration of the coloured liquid were 
found to be Z, = 75-0em and Z, = 24-4cm. When these values, together with the 
measured values of h and a, were substituted into (38), the result A = 1-85 x 10? 
was obtained. 

From the given dimensions of the cross-section of the rectangular channel, the 
ratio L = 21-1. Substituting this value into equation (36) gives A = 1-77 x 10-2, 
which is about 4° lower than the experimental value. The first term of (36), 
which corresponds to the result derived from the simple theory of §2, gives 
A x 1:66 x 10-*. This is about 10° lower than the experimental value. 

Although the contrast in the magnitude of the errors is not very great (4-10 %), 
it seems reasonable to explain the first of these in terms of experimental error, and 
to attribute the increase in the second of these to shortcomings of the simple 
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approximate theory. In particular, the measured value of h could be in error by 
about 1°%, and the inaccuracy of the cell construction may have been a contri- 
buting factor. It is also possible that the measured depth of penetration of the 
upper liquid into the Hele-Shaw cell was underestimated due to a slow progressive 
loss of coloration of the methylene blue. This effect can be produced by oxidation 
by incident light, by reaction with other substances in solution, or by bacterial 
action. A further possibility, that the denser coloured fluid ’ overshot’ its equi- 
librium position in the capillary tube, has been considered, but this seems to have 
been unlikely in the present instance, when the density contrast (p,/p) + 1-002) 
was very small. 

Figures 3, 4 and 5 (plate 1) illustrate the motion of the coloured liquid in 
the vertical rectangular channel while convection was in progress. In the early 
stages of the motion, a finger of coloured fluid passed down the middle of the 
channel, the interface near the leading edge being quite sharp (figure 3(a)). 
Before molecular diffusion had widened the coloured region appreciably, its width 
appeared to be equal to about one-half the width of the channel. Evidently the 
two liquids were behaving in a ‘quasi-immiscible’ manner, the finger of coloured 
fluid resembling those fingers observed by Saffman & Taylor (1958) in their 
experiments on the motion of an interface between immiscible liquids in a Hele- 
Shaw cell. 

Figures 3(b) and (c) show the break-up of the symmetrical finger of coloured 
fluid. While the finger was less than about 6 cm in length, the rate of descent at the 
leading edge decreased approximately exponentially from an initial rate of about 
3cm/h. However, when a depth of 6 cm had been exceeded, a fairly abrupt change 
occurred in the character of the motion, and thereafter the rate of descent was 
almost constant at about 0-09 cm/h. 

The change in character of the motion at a depth of 6em would be expected 
from the results of the stability theory. Since the maximurn depth of penetration 
was about 24cm, when the lowest mode of disturbance would be just stable, it 
follows that 6 cm corresponded to that depth of penetration at which the distur- 
bance mode with twice the wave-number became stable. When the Rayleigh 
number, based upon the average vertical density gradient of the liquid in the cell, 
decreased sufficiently for the latter mode to become damped, only a finite develop- 
ment of the lowest mode of disturbance could be self-sustaining. The subsequent 
photographs (figure 4) show the new convective motion, with the heavier coloured 
fluid moving down one side of the channel. 

The sequence of figure 4 also shows that the finite-amplitude motion of the 
liquid down the side of the channel was itself unstable, the predominant down- 
ward motion appearing first upon one side of the channel, and then upon the 
other. 

In the final sequence, figure 5, the motion leading to a stable state is shown; the 
effect of molecular diffusion across the channel became predominant, and, in 
figure 5(c), the coloured fluid had diffused right cross the channel. That stage 
indicated the completion of the convective motion. 

It is hoped that a more detailed treatment of the observed convective pheno- 
mena will be given in a later paper. 
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On the flow of a conducting fluid past a 
magnetized sphere 
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In the steady flow of an incompressible, inviscid, conducting fluid past a 
magnetized sphere, the first-order effects of the magnetic field and the con- 
ductivity are studied. Paraboloidal wakes of vorticity and magnetic intensity are 
formed, the former being half the size of the latter. The vorticity, generated by 
the non-conservative electromagnetic force, is logarithmically infinite on the 
sphere. For the case of a dipole of moment. at the centre of a sphere of radiusa, 


the drag coefficient is 144y"2 
+’) 


where and jv’ are the permeabilities of the fluid and sphere, respectively, / is the 
ratio of the representative magnetic pressure ju.M?/2a® to the free-stream dynamic 
pressure, and F,, is the magnetic Reynolds number. 


1. Introduction 


So far in the study of the effects of a magnetic field on the flow of a conducting 
fluid past a finite body, attention has been focused on problems in which the 
magnetic field is applied externally. It is of some interest to consider a case in 
which the magnetic field originates in the body itself. In all such questions the 
generation of vorticity by the action of the non-conservative electromagnetic 
force is of special concern. 

Here we consider the steady flow of an incompressible, inviscid fluid past a 
sphere of arbitrary conductivity, in which there is an arbitrary, axially sym- 
metric, magnetic distribution. For such a configuration the current lines are 
circles. The results are described for a general distribution and are evaluated 
explicitly for a dipole situated at the centre. 

Our investigation is concerned with the first-order effects of the magnetic field 
and the conductivity of the fluid. This leads to a regular perturbation in /, the 
ratio of a representative magnetic pressure to the free-stream dynamic pressure. 
A similar perturbation in the magnetic Reynolds number RF), is not uniformly 
valid; the perturbation being singular at infinity. 

The situation is similar to that in Oseen’s approximation to the viscous flow 
past a sphere: wakes (of magnetic intensity and vorticity) are formed whose 


* Now at University College London. 
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boundaries are paraboloids with linear dimensions proportional to 1/R,,. The 
vorticity wake is half the size of the magnetic wake. 

The first term in the /-expansion for the magnetic field can be obtained as a 
regular perturbation in R,, after an exponential factor (representing the non- 
uniformity at infinity) has been extracted. Once this is determined, the first term 
in the £-perturbation for the vorticity can be evaluated by integration along the 
streamlines of the potential flow past the sphere, using the fact that the vorticity 
vanishes at infinity upstream. Now the time during which a fluid element is under 
the influence of the non-conservative force, in its passage through the neighbour- 
hood of the front stagnation point, tends to infinity logarithmically with its 
closest approach to that point. Consequently, the vorticity becomes logarithmi- 
cally infinite on the sphere. 

Drag on the sphere derives from two sources, the fluid pressure and the Maxwell 
stress. For the dipole case there is no contribution from the pressure, and the 
drag coefficient is [144y’?/5(2” + where and are the permeabilities 
of the fluid and sphere, respectively.* 

The flow is completely independent of the conductivity of the sphere, provided 
it is finite. On the other hand, if it is assumed that the field was originally frozen 
into the sphere, the results are different. This case is discussed in the last section. 


2. Basic equations 


The equations governing the steady motion of an incompressible, inviscid, 
electrically conducting fluid of constant properties are: 


(a) divv = 9, (b) curlvxv = —grad [(p/p,)) + 407]+ (u/p,) curl H x H, 
(c) curlH = o(E+yvxH), (d) curlE=0, (e) divH=0. (1) 
For an axially symmetric motion in which the velocity v and magnetic intensity 
H lie in the meridian plane and are independent of the azimuthal angle, the 
conduction equation (1c) shows that the electric field E is perpendicular to this 
plane and also independent of the azimuthal angle. It then follows from (1d) that 
E = 0, if it is to be finite at the axis. When v, r and H are now made dimensionless 
by referring them to the velocity at infinity U, the radius of the sphere a, and a 
representative magnetic intensity h, respectively, the equations reduce to 


(a) divv=0, (b)curlvxv = —gradP+/curlH x H, as 
(c) curlH = Ry, vxH, (d) divH = 0, (?) 
bh? 
where P=p+ = Ry, = Uano, 


and p is now the pressure divided by p, U?. 

When there is no magnetic field, / is zero and (2b) shows that P is constant 
along streamlines. From the assumed uniform conditions at infinity it then 
follows that P is constant throughout and hence [using (26) again] that curl v = 0. 
In conjunction with (2a) this yields the potential solution 

v= v= |(1- 5) sino, oF (3) 


~ 


* The total drag is more easily obtained from the Joule dissipation (Chopra & Singer 1958). 
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satisfying the condition of zero normal velocity at the sphere r = 1 and tending 
to the given uniform stream (cos9, —sin@, 0) at infinity. Here we use spherical 
polar co-ordinates with 9 = 0 pointing downstream (see figure 1). The corre- 
sponding pressure field p = p, is then determined by the constancy of P. 

For weak magnetic fields (small h) we expand in powers of £ 


=VotAv,+.... H=H,+fH,+.... 


Note that H, is not zero. To obtain the magnetic field, H must be multiplied by 
h and it is the latter which tends to zero in the limit. According to (2c) and (2d), 


curlH, = Hy, div H, = 0, (4) 


= Const. 


FicureE 1. Integration along the potential flow streamlines. 


where V, is given by (3). The remaining equations (2a) and (25) give 
curly, x V, = —grad P,+curlH,x Hb), div v, = 0, (5) 
and these determine v, and p, once H, is known; here 
= Pi+Vo-V1 


and we have used the fact that curl v, = 0. 

We shall restrict the discussion to Ho, v,, and p,. It is easily seen, however, 
that at each stage in the approximation the terms in v and p are determined 
before the term in H. 


3. Determination of H, 


The subscripts 0 on H, and 1 on v, and p, will now be dropped. In axially 
symmetric flow the second of equations (4) allows us to write for H = (H,, Hg, 0) 
1 1 0A (6) 


00’ rsinO or’ 


where A is a function of r and 4 alone. On introducing this into the single non- 
zero component of the first equation, we find that A must satisfy 


( : ssa) = Ry| (1-53) sind. (7) 


r= or or 273 00 
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For small values of R,, an appropriate solution of this equation may be found 
by a perturbation expansion. The perturbation is not regular at infinity, however, 
since the right-hand side of (7) vanishes more slowly than the left, as r > oo, when 
A is algebraic. In fact, if small terms in r on the right-hand side are neglected, the 
equation becomes* 


of which simple solutions, for arbitrary R,,, are 
A = exp cos sin? (cos 4). (8) 


Here n is a positive integer, K is a Bessel function, and P the Legendre poly- 
nomial. Apart from the exponential, A is a product solution, chosen so as to be 
regular on the axis and finite at infinity. For large values of r, 


A ~ exp — cos @)]sin? 6 (cos 4), 


which gives an exponential decrease for all 6 + 0. On the other hand, for r fixed, 
tends to a multiple of P’,(cos 0)/r” as Rj, > 0; that is, a solution of 
(7) with Ry, = 0(curlH = 0). 

Thus the magnetic field is swept into a wake behind the sphere, whose boundary 
is the paraboloid of revolution r(1 — cos @) = ¢/R,, (c a suitably chosen constant). 
As R,, > 0 the wake broadens so as to fill the whole of space in the limit. 

It is now clear how the perturbation must be carried out. In equation (7) we set 


A (p = Ryr), 
and obtain for « the equation 


op? cp p® 


1 
= — Ri, |— (1 (1 +3080) (9) 


op’ 2p* 
With R,, = 0 this has the product solutions 
a = ..(4p)sin? (cos 0), (10) 


in accordance with (8). We now take a suitable combination of these solutions (so 
as to satisfy the boundary conditions at the sphere) and use it to generate the 
solution of (9) by successive approximation, ensuring at each stage that the 
boundary conditions are satisfied (to the appropriate order in R,,) by adding a 
suitable combination of (10). 

The factor #3, on the right-hand side of (9) does not imply that the perturbation 
is O(R4,). On replacing p by Ry,r, it becomes 2,,; and at any finite r the approxi- 
mations proceed in powers of R,,. However, there is an immediate exception due 


* In Oseen’s approximation to the viscous flow past a sphere, this is the equation 
satisfied by the vorticity, which is hardly surprising since the same physical processes 
are involved and the same type of approximation has been made [v, replaced by the 
velocity at infinity in (4)]. 
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to the fact that the disturbance part of v, is a dipole field: the first two terms on the 
right in (9) cancel when the leading term in the first approximation is substituted. 


The correct combination of (10) to take when there is a dipole at the centre of 
the sphere is 


a =KR,,|—+ sin? 0+ ARS, sin?@ cos0, 

containing the solutions n = land n = 2.* Herex and A are constants which must 
be determined in conjunction with the field inside the sphere. Since the first 
perturbation, found by substituting this expression into the right-hand side of (9) 


and integrating, provides corrections O(R%,), we may neglect the uncorrected 
terms at this stage and set 


A = exp[—4R,,7r(1—cos4)]{(k/r) sin? 6 + Ry,[4« + (A/r?) cos A] sin? 6}. (11) 


~ 


FicureE 2. Distortion of the magnetic field of a dipole, for Ry = 0-1 and yw’ = uw. Broken 
lines show undisturbed field. Inside the sphere (r = 1) the fields are indistinguishable. 


Inside the sphere we again have E = 0, so that curlH = o’E = 0, whatever 
the conductivity o’ of the sphere. Hence we take 
A = (1/r) sin? 6 + (kr? + cos @) sin? 4, (12) 


the first term being due to the dipole and the rest representing the (regular) 
irrotational disturbance field (k and l are constants to be determined). We can now 
identify h as M/a®, where M is the moment of the dipole. 

The validity of (11) and (12) now follows from the fact that «, A, k and 1 can be 
chosen so that the boundary conditions on the sphere are satisfied to order R,,. 
Thus the continuity of the tangential component of H requires 


1—2k=x, 3l= —2A, 
and the continuity of the normal component of magnetic induction requires 


wl=p(hx+A), 


* K,,.4(z) is a polynomial of degree n in 1/z multiplied by e-#/zx}. 
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where yu’ is the permeability of the sphere. These are the only conditions; and they 
lead to the values 
We shall not take the approximation any further than (11) and (12), though it 
may easily be done. From (6) we find for the field outside the sphere 


_ 


Dee 
H, = exp[ (1+3R,,r)cosO+ Ry + 3A) cos? — + 
(13) 


H, = exp[—4R,r(1 —cos0)}!“ (14 $Ryr) sin — 2A)sin0 cos), | 


\r3 


to the same degree of approximation (see figure 2). 


4. Determination and behaviour of the vorticity 


The vorticity vector curl v is perpendicular to the meridian plane, and hence 
has only one component, w. According to (5), it satisfies the equations 


+ sin@ = tk, (14a) 
1oP 

o(1-5) cos @ = ag t (146) 


where F = (F,, Fj, 0) = curlH x H = R,,(v) x H) x H. From (3) and (13) we find 
F, = exp[—R,,r(1—cos 0)] { — (3«?R,,/r®) (1+ sin? 6 cos 
+ sin? O[(3Kr> — L1Ar? — 3xr2 + 4A) cos? + (273 + 1) (fxr? + JA)T}, 

(15a) 

F, = exp[— R,r(1 — cos 6)] {(6x?R,,/r°)(1 + sin cos? 

+ (KR5,/r!) sin cos O[(3xr® + 19Ar3 + — A) cos? 6 — (5r? + 1) +. A)]}, 

(155) 

so that on eliminating P from equations (14), we obtain 


1 0b 
cos (1 + sa) = f(r,0; Ry), (16) 
where © = w/(rsin@) and 


rsin@| cer 


= (1+ Ry,7r) (7 cos? 6 +1) 


2 
cos O[( — — 70Ar3 — 12 Kr? + 7A) cos? + (9Kr? + + + 


Note that the argument of the exponential has been doubled, which means that 
the size of the vorticity wake is only half that of the magnetic field. 
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The left-hand side of (16) represents a differentiation along the streamlines of 
the potential flow (3). Hence it is convenient to introduce the potential function 
and stream function of this flow 


1 


When these are taken as new independent variables in place of r and @, the 
@-equation becomes 


(1-1) sinta 0:R 


0 


Thus we find that at any point P 


op r f(r,0; Ry,) 
¥ 1)? cos? 6 + (r? +4) 2 sin? (17) 


where the integrand is to be expressed in terms of ¢, and yo, and the integration 
taken with yj, fixed. The choice for the lower limit ensures that on each streamline 
w > 0 as r-> oo upstream. This is the only boundary condition needed to deter- 
mine the vorticity. 

Of special interest is the behaviour of w at the sphere. Note that the integrand 
of 6 becomes infinite at A (figure 1), through which passes the dividing streamline 
connecting —0o to any point on the sphere. Hence we consider the integration 
over a short segment P, P, of a neighbouring streamline (yf, small), and let y» tend 
to zero (keeping the limits d, = ¢,, ¢2 fixed). Now, for wy) = const., 


dr 
= [(73— cos? 3.4.1)2 gin? 
= [(r3— 1)? cos? + (r3 + 4)? sin — 1) e080” 


so that the dominant part of the integral is 


I= “| dr = (1,7; log 
the remainder wailies to a finite limit as 7, > 0. Since the integrals taken over 
the ranges (— 0, ,) and (4, dp) also tend to finite limits, this is the leading part 
of the complete integral (17)+. Finally, we note that 7,—1 may be replaced by 
rp—linJ,since P,and P lie on the same streamline. It follows that, if terms which 
tend to finite limits are omitted, 


w = —}sindpf (1,7; Ry) log (rp—1) 
= + 2A) Ry, + O(R5,)] sin Op log (rp— 1), (18) 


when P is close to the sphere. 

The vorticity is logarithmically infinite only on the sphere. For a point P close 
to the axis downstream, a term log (r, — 1) still arises from integration near A (and 
a similar term from integration near the rear stagnation point B), but this time 
it is replaced by log (sin? ,,), a quantity which, when multiplied by sin 4, tends 
to zero with 7,. It is easily seen that w is also zero on the axis upstream. 


* This is the square root of what is usually taken to be the stream function. 
+ For more details, see Appendix 1. 
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It is also of interest to see how the vorticity dies off at large distances. More 
precisely, we shall assume that ¢, = R,,rp(1—cos@,) is large: in particular, this 
ensures that r is large at every upstream point on the streamline through P, so 
that the integrand in (17) may be replaced by 

RS 
(cos? 0 +7 cos? 0 — 3 cos 0 +1) exp [— Ry, — cos 


> 
‘ 


while, effectively, =reos0, wy =rsind. 
If now we change the integration variable to 


r(1—cos@) 
rp(1—cos 


the integral may be written as 


~ 


=“ 


e§P"g(u, du, 
P 0 


where g is a rational function of u with 
g(0, = (cos? Op +7 cos? Op — 3 cos Op + 1) 


and poles at wu = —1+icot(40,). It now follows from Watson’s lemma [see 
Copson 1935, p. 218] on asymptotic expansions that, for each fixed 6p, 
-2 


— Sindp (cos? 6p + 7 cos? 0p — 3 cos @p + 1) 
P 


x exp[—Ryrp(1—cos p)] +0(z-)| 


when &>p is large. 


5. Expansion of w 

The velocity field corresponding to can be written down by using the known 
velocity field of a circular vortex in the presence of a fixed sphere (Yeh, Martinek 
& Ludford 1955). The resulting formula is too complicated for our purposes, and 
accordingly we resort to a series expansion for w, which, although it can in 
principle be derived from (17), will in fact be obtained directly from (16). It 
turns out that only certain (Fourier) components of the velocity field can be so 
determined, but fortunately the only component needed for the computation of 
the drag is one which can. 

The stream function yy corresponding to satisfies the equation 


ey sindd (1 oY 


and in terms of it the velocity components are 
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Since the product solutions of (19) for o = 0 which are regular on the axis have 
the form 


(4, sin? (cos) (n a positive integer), (21) 
it is clear that we should expand in the form 


w= w,,(r) sin OP, (cos 6). (22) 


n=1 


Correspondingly the right-hand side of (16) is written 


f(r,0; Ry) = Q,(r) 9). (23) 
n=1 
where the 2, are determined by suitable expansions of F, and £; (see next section). 
When the series (22) is substituted into the left-hand side of (16) and the result 
brought into the same form as the right-hand side (23), we find the recurrence 
relation 


n+2 [/ 1 = 
n+3 (r- + (» +2 +— On 


on equating coefficients of P’,(cos 4). Its integrated form is 


n+2 (r3—1)hm+3) Pr —] n—l1T/ 


The lower limit of integration is determined by the requirement that w, should 
not become algebraically infinite on r = 1. 
For n = 1 we obtain 


ar 


= (a — 1) QO, (a) da, (24) 
and then 4, W¢, s, .... follow nat On the other hand, for n = 2 the 
relation involves both and w;, and w, must be chosen before ... Can 
be determined. This must be done in such a way that the complete series (22 
tends to zero at infinitely large distances upstream, and, presumably, , is fixed 
by this condition. In any event, in order to obtain the drag on the sphere it is 
sufficient to know w,. 
The y,,(r) in the expansion 


y= y,(r)sin? 0) 


of the stream function is determined, according to (19) and (22), by the equation 
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The solutions which vanish on r = 1 are given by 


l r Eats 1 
__ 
and we choose A, = mal [ gn-t dé, (26) 


so that y,, will remain finite at infinity. 
Once the y,, have been determined, the velocity field follows from (20). In 
particular, 
= v,(r)sinOP, (cos), v,(r) = — Wnl(r)/r. 
n=1 
The only detailed piece of information we need concerning this series is the value 
of v, on r = 1 (see §6). From (25) and (26) we find 


= —Y,(1) = —5A, 


~ 


5 [2 


Q,,(a) dx. (27) 
9], 


6. The drag 


The force exerted by the fluid on the sphere has two parts, one due to the fluid 
pressure and the other to the Maxwell stress. We consider these in turn. 
Assume that the F, appearing in equation (14) has been expanded in the form 


a,,(r) sin (cos ). (28) 
n=1 


Then on the sphere we find by integration of (145) that 
P = p+vy.v = — ¥4,(1)P, (cos), 
n=0 


where a,(1) is an undetermined constant of integration. Hence the total drag 
force due to pressure is p,U*a?D,, where 


= pcos). 2m7sin 


Pp 0 
7 
= | |- 3 sind a,(1) P,,(cos 0) sin cos @ dé 
0 n=1 n=0 
= —§v(1) + 
That the remaining terms in the series integrate to zero follows from the ortho- 


gonality property of the Legendre polynomials. 
To complete the calculation we assume that F, has been expanded in the form 


b,,(r) P,,(cos 8). (29) 
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Then the Q,, appearing in the expansion of f, see (23), may be written 


and the integral (27) becomes 


Thus we finally have D, = imp |” da. (30) 
1 


So far the treatment has been quite general. For the present case of a dipole 
field we find 
D, = (31) 


Details of the calculation can be found in Appendix 2. 
There remains the drag due to the Maxwell stress, wH; H; — $H?6,;, which on 
a surface element of the sphere has components 
wh?{\(H2— H3), H, Hy, 0}. 


The total contribution is a force p, U?a?D,, where 
Dy, = cos 0 — H, H,sin @]sin 6 dd. (32) 
Jo 


In the present case the square bracket is 

4x2(4 cos? — 5 sin? 7) cos 0 + KRy,[(4« +A) (3 cos? — 1)? — 6A sin? 6 cos? 6] + O(R5,) 
[see (13)], and, as we expect, the first term makes no total contribution to the drag. 
The remainder give Dy = + (33) 


Note that the formulas (30) and (32) enable one to calculate the drag directly 
from the magnetic field. Also note that (33) is correct to quadratic terms in / 
and R,, since if either / or Ry, is zero there can be no drag. 


7. The perfectly conducting sphere 

The preceding analysis applies to a sphere of arbitrary conductivity, and hence 
to one of arbitrarily large conductivity. However, if the conductivity is infinite, 
so that the magnetic field has been frozen into the sphere, the results are different. 
We consider in detail the dipole case again. 

The field outside is again determined by the solutions n = 1, 2 in (10); neglecting 
uncorrected terms as before, we find 

1 


A = sin? 0+ 0) sin?) (34) 


On the sphere A = (y’/) sin? 0, correct to O(R,,), so that the normal component 
of magnetic induction [see (6)] has the same value as that for the dipole field 
frozen into the sphere. There is no requirement on the tangential component of 
the magnetic field; any discontinuity in it corresponds to surface currents on 
the sphere. 
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The results are now obtained without further calculation by noting that (34) 
can be derived from (11) by formally setting « = —2A = y’/u. Thus, according to 
(18), the vorticity again becomes logarithmically infinite at the sphere: 


when P is close to the sphere. As before, there is no contribution to the drag from 
the pressure; the total drag (due to the Maxwell stress) is py) UZa*D,,, where [see 
30) equation (33)] 
8m 
Dy = — (“) 
ole 5 \u 
31) From this we see that when a sphere of high conductivity is considered to be 


perfectly conducting, the computed drag is in error by a factor 


+h 
For example, soft iron has a w’ of the order of 300 for weak fields, and the factor 
is about 10-4. In equation (12), k = —1 andl = 0 for large y’/u, so that the lines 
of induction tend to become compressed into the sphere. On the other hand, the 
32) frozen-in field approximation would be appropriate for materials used in per- 
| manent magnets; it would also hold for a uniformly magnetized sphere. 


We are indebted to S. Goldstein for suggesting this problem and for his interest 
i) at every stage. The research by one of the authors (G.S.S. L.) was supported in 
ag. part by the United States Air Force under Contract no. AF 49(638)-154, moni- 
tored by the AF Office of Scientific Research of the Air Research and Develop- 


33) ment Command, and in part by the Office of Ordnance Research, U.S. Army, 
tly under Contract DA-36-034-ORD-1486; and that by the other (J. D. M.) by the 
n p United States Air Force under Contract AF 19(604)-4545 with the Geophysics 
Research Directorate of the Air Force Cambridge Research Center, Air Research 
and Development Command. 
nce Appendix 1 
ite, In view of the surprising nature of this result, we feel a brief discussion of the 
mt. analysis should be given here. 
; For ¢, and ¢, fixed, it is clear that the integral (17), taken over the partial 
ing ranges (—0o, ¢,) and ¢,), tends to finite limits as — 0. Hence we must 
how that 
show tha Pf (7,0; Ry) 
34) ~ Jy, cos? + + 4)? sin?6] *° 
ail also tends to a finite limit. 
eld Now, in a sufficiently small neighbourhood, ¢, < $ < $2, 0 < Wo < €, of the 
of front stagnation point, we can find constants A and B such that 
on | r2f(r,0; Ry) f(1,m; Ry) 
= —1)+B-,—*— sin 0. 
|(r2+r-+ 1) (-—3) r(r+r+1) 


| 
f 
2 
= s( Ry, {1+ O(R5,)]sin log (rp—1) 
‘ 
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Also, for yp fixed, 
dr rdé 


(r3— 1)? cos? + (r3+4)?sin?0 (r+. 4) sind” 


<| K cos dé 


= A(r,—r,)— B(sin —sin 


Hence |J 


which clearly tends to a finite limit as yy —> 0. 


Appendix 2 
From equations (28) and (29) we have 


0 


= x sin + 2F, cos sin dé, 
0 


where, according to equations (15), 


F, sin 0 + cos 


= exp cos O[(3Kr? — A) cos? — (4x7? +A)]. (35) 


Note that the leading terms in F, and F, cancel, this being a peculiarity of the 
dipole field. 

Now, in evaluating D,, from (30), we may set Rj, = 0 in this last exponential, 
since this leads only to an error O(R4,). But then (35) is antisymmetric about 
= 4 iv 
) = 4m and hence gives a,(r)-+By(r) = 0. 


From this follows (31). 
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The Stokes flow problem for a class of axially 
symmetric bodies 


By L. E. PAYNE 
Institute for Fluid Dynamics and Applied Mathematics, University of Maryland 


AND W. H. PELL 
National Bureau of Standards, Washington 25, D.C. 


(Received 6 July 1959) 


The Stokes flow problem is concerned with fluid motion about an obstacle when 
the motion is such that inertial effects can be neglected. This problem is con- 
sidered here for the case in which the obstacle (or configuration of obstacles) has 
an axis of symmetry, and the flow at distant points is uniform and parallel to 
this axis. The differential equation for the stream function yy then assumes the 
form L?.,y = 0, where L_, is the operator which occurs in axially symmetric 
flows of the incompressible ideal fluid. This is a particular case of the fundamental 
operator of A. Weinstein’s generalized axially symmetric potential theory. Using 
the results of this theory and theorems regarding representations of the solutions 
of repeated operator equations, the authors (1) give a general expression for the 
drag of an axially symmetric configuration in Stokes flow, and (2) indicate a 
procedure for the determination of the stream function. The stream function is 
found for the particular case of the lens-shaped body. 

Explicit calculation of the drag is difficult for the general lens, without recourse 
to numerical procedures, but is relatively easy in the case of the hemispherical 
cup. As examples illustrative of their procedures, the authors briefly consider 
three Stokes flow problems whose solutions have been given previously. 


1. Introduction 

The determination of the isothermal flow of an incompressible, viscous fluid 
about an impermeable body immersed therein requires the solution of the Navier— 
Stokes equations and the equation of continuity subject to the condition that the 
velocity of flow coincide with that of the external boundary of the body at each 
of its points. The non-linearity of the Navier-Stokes equations renders the 
solution of this problem extremely difficult, and only a relatively small number of 
exact solutions of rather specialized character are known (Dryden et al. 1932; 
Lamb 1932). In many instances, however, the flow is such that plausible 
simplifying assumptions regarding its character (in certain regions. at least) 
can be made which result in a less refractory mathematical problem. A large part 
of the theory of viscous flows consists in the discussion of problems obtained 
in this manner. 

The oldest problem of this type is the so-called ‘Stokes flow’ problem. It is 
defined by the assumption that inertial effects are negligible in comparison with 
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those of viscosity, or, more precisely, that the Reynolds number R of the flow is 
very small. This situation obtains when either the characteristic flow velocity or 
body dimension (or both) appearing in FR is suitably small or the kinematic 
viscosity is large. Sir George Stokes (1850), in the course of treating the steady 
motion of a sphere in a viscous liquid, appears to have been the first to omit the 
inertial terms in the equations of motion. Oberbeck (1876) solved the Stokes flow 
problem for an ellipsoid immersed in a uniform flow, and also treated the special 
cases of a circular disk broadside and edgewise to the undisturbed flow. More 
recently, Ray (1936), by a more direct method, has obtained the solution of this 
circular disk problem in a different form. His solution yields the same drag as was 
obtained by Oberbeck. The problem of a sphere moving slowly near a plane wall 
was solved by Lorentz (1897), while that of a sphere falling along the axis of a 
vertical tube has been discussed by Ladenburg (1907). 

A number of two-dimensional flows of Stokes type have also been treated. 
Berry & Swain (1923) have considered the elliptic cylinder immersed in an 
infinite body of fluid, as well as the special cases of the circular cylinder and the 
infinitely thin flat plate athwart and edgewise to the flow. The speed of flow 
becomes logarithmically infinite at points infinitely far from any of these bodies. 
Dean has devoted a series of papers (1944; others are referenced there) to the 
Stokes flow past protuberances of various shapes in otherwise flat walls, the fluid 
filling the half-space lying on that side of the wall which bears the protuberance. 
Here again the speed becomes infinite at an infinite distance from the wall. Dean 
(1936) has also considered the Stokes flow problem for a fluid which emerges 
uniformly from an aperture in a plane wall into a semi-infinite body of fluid. In 
this case the speed of flow vanishes at infinity. Green (1944) has found a Stokes 
flow for the region bounded by converging plane walls and for the region lying 
between two branches of a hyperbola. Again the speed vanishes at infinity. 
Further discussion and bibliography on Stokes flows may be found in Dryden 
et al. (1932) and Lamb (1932). 

It seems reasonable to doubt the uniqueness of the Berry—Swain solutions, 
since it appears not unlikely that there may exist a Stokes flow about the cylinder 
which is uniform at infinity, in view of the Stokes and Oberbeck results for the 
sphere and ellipsoid. It has been generally accepted, however, since the time of 
Stokes, that such is not the case, Stokes himself giving a physical argument for 
this apparent anomaly (1850). The question has recently been settled once for all 
by Finn & Noll (1957), who have given a careful proof of the fact that the only 
Stokes flow uniform at infinity about a cylinder whose boundary is composed of 
a finite number of piecewise smooth non-intersecting simple closed curves lying 
in the finite plane is the state of rest. The boundaries in the other two-dimen- 
sional problems mentioned above do not lie entirely in the finite plane, and 
so the Finn—Noll theorem does not apply. 

In the case of three-dimensional flow, the prescription of a uniform velocity at 
infinity in a Stokes flow gives a well-set problem, for Finn & Noll (1957) have 
demonstrated the uniqueness of such a flow about a body whose surface is com- 
posed of a finite number of piecewise smooth non-intersecting simple closed 
surfaces. See also Odqvist (1930). 
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In this paper we propose to consider the solution of the Stokes flow problem for 
axially symmetric bodies with the aid of the generalized axially symmetric 
potential theory initiated and developed by Weinstein (1948, 1955) and his 
co-workers (Huber 1953; Hyman 1954; Payne 1958). After the formulation of the 
problem and a general discussion of the techniques employed in its solution, we 
restrict our discussion to simply connected regions (in a meridional plane). We 
first obtain a general expression for the drag of the simply connected axially 
symmetric body, and then solve the Stokes flow problem for the general lens- 
shaped body. Some interesting special cases of this body are discussed. As 
examples of the method, some previously solved problems are discussed briefly. 
Finally, a table of the drag of various bodies in Stokes flow is given. 


2. Statement of the Stokes flow problem for axially symmetric bodies 

Consider a collection of n bodies which are individually axially symmetric and 
which are so arranged that the same is true of the aggregate of bodies. Let these 
be immersed at rest in a uniform flow of a viscous, incompressible fluid which 
fills three-dimensional space, and let the axis of symmetry be parallel to the 
direction of the uniform flow (referred to as the free stream direction). We refer the 
flow to cylindrical co-ordinates (x, r,). The a-axis is chosen to lie along the axis 


C; 
(x,r) 
NN C 
O x 


FicuRE 1. The general configuration. 


of symmetry of the body with its positive direction the same as that of the free 
stream, r is radial distance from this axis, and @ is an azimuthal angle defining 
meridional planes through the x-axis. In view of the body-flow—co-ordinate 
configuration the co-ordinate @ will play no role in our analysis, and we may 
restrict attention to any single meridional half-plane, as shown in figure 1, so that 
the range of co-ordinates is 27. A meri- 
dional section of the bodies consists of regions bounded by two types of curves: 
(i) closed contours C;, i = 1, 2,...,1, each of which lies entirely above the x-axis, 
and (ii) ares C;, 7 = 1,2,...,k (k+l =n), which have termini A, and 6,, but no 
other points, lying on the x-axis (figure 1). The region of flow D is that portion of 
the meridional half-plane r > 0 which is exterior to the regions enclosed by the 
Cand the C; + A, B;. The C,; and C; are assumed to be made up of a finite number 
of smooth ares, joined end to end. For brevity, we shall refer to a meridional 
34-2 


iS 
r 
: 
ul 
IS 
iS | 
a 
n 
2, 
n 
> 
n 
Ve 
n 
4 
3, 
| 
of 
| fe 
\ 
of | 
d 
d 
| 


532 L. E. Payne and W. H. Pell 


section of a body as ‘the body’, and the bounding curve C, or C; as its profile. The 
complete set of profiles constitutes the boundary of the set of bodies and is 
denoted by C. The complete boundary of the flow region, however, consists of 
C together with those segments of the x-axis exterior to the set of A;B,, 
j= 1,2,...,k. The velocity of flow at a generic point (x,r) of D we denote by 
u(z,r) = (u,(x,7r), u,(x,r)), and the undisturbed (free stream) velocity by 
u,, = (U,0), where U > 0 is a constant. We assume that 


lim u(z,r) = u,, (2.1) 
where p = (a2+r?)h, (2.2) 
Since the fluid is incompressible and the flow steady, the continuity equation 
becomes simply 


The axial symmetry of the flow permits the introduction of a stream function y. 
Accordingly, we write 
lov 

4 = (2.3) 


and the continuity equation is then of necessity satisfied. 

The vorticity vector ¢ = curlu has at each point of the meridional plane the 
form (0,0, ¢) and it is thus sufficient to deal only with the scalar ¢. It is seen at 
once (Milne-Thomson 1950, p. 494) that 


Cu, CU, 


(2.4) 
and if (2.3) are inserted in this we obtain 

1 9 

Liyy, (2.5) 


q2 

10 

where L_y = (2.6) 
ox? Or? 


The reason for the departure from the usual notation (Milne-Thomson 1950, 
p. 494) for this operator will be clear later. 

By a well-known procedure (Milne-Thomson 1950, pp. 507-9) the Navier-— 
Stokes equations can be converted into the following equation for ¢ in D: 


0 
(Cu,) += (Lug) += = 0, (2.7) 


where v = s/p is the kinematic viscosity. Substitution from (2.3) and (2.5) then 
permits us to write this in the form 


oly, 


2. 
=0. (2.8) 


The Stokes assumption asserts that the first term of this is negligible in com- 
parison with the second, and we thus obtain for axially symmetric Stokes flow 
the equation Day =0 (2.9) 
to be satisfied in D. 


c 


so 


in 


u 
t 
t] 
St 
(a 
in 
3. 
ar 
sa 
CO 
eq 
de 
a 
al 
th 
| = Mc 
cal 
It 
th 
de] 
sat 
for 


Stokes flow for axially symmetric bodies 533 


The kinematic condition of vanishing normal and tangential components of 
u at the boundary of an impermeable body at rest in a viscous fluid leads to the 
conditions 


y = k, a constant, (2.10a) 
(2.105) 
On 


to be satisfied on each C; and C}. We take n to be the unit normal on C directed into 
the fluid. In general, k will be different on different C, although it will have the 
same value on all C;,. 

Finally, it is easy to see from (2.1) and (2.3) that y% must satisfy the condition 
lim y(x,r) = 3r°U (2.11) 

(apart from an arbitrary additive constant which we have taken to be zero). 

In summary: the solution y(x,r) of (2.9) subject to the conditions (2.10) and 
(2.11) solves the problem of Stokes flow for an axially symmetric body immersed 
in an infinite fluid whose velocity at infinity is parallel to the axis of the body. 


3. Representation of the solution of the flow equation 


We follow the usual procedure in solving boundary value problems: we obtain 
solutions of the differential equation (2.9) which seem suitable for our purpose, 
and construct therefrom by the well-known linear operations one which also 
satisfies the boundary conditions (2.10)-(2.11). The writers believe the problem 
considered here to have some novelty, however, inasmuch as the differential 
equation (2.9) is one which has received little attention and also because the 
solution given here constitutes an application of powerful methods recently 
developed by Weinstein (1948, 1955), Payne (1958), and other workers (as 
indicated below) in generalized axially symmetric potential theory. 

Let y*(x,r) denote any solution of 


kov 
Or? ror 4) 
(—o« < k < w) in D which is regular on the x-axis outside of C. We refer to y* as 
a generalized axially symmetric potential function. The behaviour of such 
functions in the neighbourhood of the z-axis has been discussed for all real k by 
Hyman (1954) and Huber (1953). It can be checked by direct substitution that 


the following functions are solutions of equation (2.9) for the stream function: 
(a) (b)  (c) p?r?y?, (d) (e) rtp. (3.2) 


Moreover, Payne (1958) has shown that in certain regions any solution of (2.9) 
can be represented as a linear combination of any two of the expressions (3.2). 
It may be surmised (and this will be verified in the examples we shall consider) 
that the optimal combination to be chosen as a solution for a specific problem will 
depend on the geometry of the boundary C on which (2.10) and (2.11) must be 
satisfied. A combination suitable for one problem may be completely intractable 
for another. In fact, we choose quite different combinations in discussing the flow 
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about a spheroid from those used in the case of the lens-shaped region or separated 
spheres. 

In the remainder of this paper we shall restrict our attention to the case in 
which only profiles of the type C; occur. The connected curve formed by the 
C;,j = 1,2, ...,k, and those segments of the x-axis exterior to the A, B; (figure 2) 


FicureE 2. Configuration for simply connected flow field. 


then form a streamline of the flow, and, in fact, constitute the boundary of D. 
Thus D is simply connected for this configuration. Since y = 0 on the z-axis at 
infinity, it must remain so along the entire streamline above, i.e. on C and the 
z-axis exterior to the bodies. Condition (2.10) holds on each C; and thus on C. 
Conditions (2.10) revised for the case under consideration are thus 


0, (3.3a) 
(3.30) 
on 


on C. 

If bodies of type C; occur, then D is multiply connected, and (3.3a) no longer 
holds. 

It is convenient to write the stream function yy in the form 


= (3.4) 


and to formulate the problem in terms of 1, rather than y. It is clear that y, 
satisfies the same differential equation (2.9) as yy. In addition, it must give rise to 
a vanishing velocity at infinity, and fulfil the conditions 


yy, = 40Pr, (3.5) 
(3.55) 
cn on 
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4, Drag on the axially symmetric body 


The drag on an axially symmetric body at rest in the flow, expressed in 
cylindricai co-ordinates, can be shown to be 


mu 2nu 


(see Milne-Thomson 1950, §19.21). 
We next show that this integral for the drag may be replaced by the expression 


To this end we use the identity 


which can be obtained from the divergence theorem by a procedure resembling 
that used in deriving the symmetric form of Green’s theorem (see Weinstein 
(1948), equation (20)). The functions wu and v, the region D* and its boundary C* 
are subject to the hypotheses of that theorem, and n is a unit normal to C* 
exterior to D*. 

We nowset u = L_,y,,v = yy, and apply (4.3) to the region D* whose complete 
boundary C* consists of (i) the ares C,, (ii) a semicircle [ defined by p = R, r > 0, 
where R is arbitrary save that all profiles C; must lie within p < &, and (iii) those 
segments of the x-axis which lie between the termini of I and are exterior to the 
bodies. Since L*., yy, = 0, we see at once that we then have 


en 1 l 
| Par de = ds. (4.4) 


Since yy, is a solution of (2.9) it follows from the significance of the notation that 
L_, ys, = Ww. We now make the plausible assumption that at all points of C on 
the a-axis the velocity gradients are bounded. It then follows from equation 
(2.5) and the results of Huber (1953) and Hyman (1954) that L_, is O(7r?) as 
r-> 0 along the x-axis outside the body. Moreover, it is apparent from either the 
decomposition formula of Weinstein (1955) or those of Payne (1958) (which are 
valid in the neighbourhood of 7 = 0) and the results of Huber and Hyman that 
if yy, vanishes at 00 on the 2-axis and is analytic in a neighbourhood of a segment 
of the z-axis it must vanish as r? there. From these remarks it follows that the 
integral on the right-hand side of (4.4) arising from that portion of C* which lies 
along the x-axis vanishes. We thus have 


(L(y) P drdx = - Lar) ds 
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which may be rewritten as 


| 
D* on on 

cer on 41 Ss. ( ) 


Since L_,(r?) = 0 and y, is a solution of (2.9) we may use (4.3) to show that the 
integral over C* above vanishes, and we obtain the important result that 


From the decompositions (3.2) given by Payne it follows that we may write 


= PUY + YP) (4.8) 


in any region sufficiently far removed from the body. Since y! and y* may, 
however, be regarded as axially symmetric harmonic functions in spaces of 3 and 
5 dimensions, respectively, (4.8) yields the representation 


(K a constant) in the neighbourhood of p = 0. Further, we may conclude from 


Wi 


this that | 


) (4.10) 
in the neighbourhood of p = ». 

With this information we return to (4.7) and let & — oo in that equation. This 
yields, since D* > Das R > 


Jo R 0 


or P= 8k. (4.11) 


In view of (4.9) K = lim pvw,/r?, and we have therefore established (4.2). We now 
P¥ 1 
turn to the consideration of a specific flow configuration. 


5. The flow about a lens-shaped body 

In order to compute the flow about a lens-shaped body we introduce peripolar 
co-ordinates (Hobson 1931, pp. 422, 451) (£, 7) in the plane of the (x, r) co-ordi- 
nates through the transformation 


é bsinh (5.1) 


~ cosh — cos &’ cosh 7 — cos 


In effect, this is a dipolar system of co-ordinates in which the family of co-axial 
circles nest about (0, +b). These circles are the curves 7 = const., while the curves 
£ = const. are ares of circles joining (0, + 5) (i.e. ares of the orthogonal trajectories 
of the circles 7 = const.). As before, we consider only the half-plane r > 0. 


: | a 

(4.7) 

anc 
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The region bounded by the surface of revolution obtained by revolving two 
ares = £,, and £ = , about the x-axis is referred to as a lens-shaped region 
(figure 3). A lens-shaped body is thus one which occupies a lens-shaped region. The 
interior of the region is defined by 

0<§<£<&<29 (©<9<a@) (5.2) 


and the exterior is then given by 


(0<79 < (5.3) 
Accordingly, our problem is to find a solution of (2.9) which satisfies the 
conditions 0, (5.4a) 
0, (5.45) 
r 
(0,6) 
(7 


(0,—6) 


FicureE 3. The lens-shaped body. 


Let & be a value of € such that 
(5.5) 


We choose for yy a representation of the form 


(s—1 
—— —(s—t)t | F(a, £) K®(s)da!, 5.6 

where F(a, £) = cos A(x) cosh a + B(x) sinh ag] 
+ sin £[C(x) cosh a€ + D(a) sinh a€]. (5.7) 


The functions A(x), ..., D(x) are to be determined by the boundary conditions, 
s=coshy, t= cosé. (5.8) 
After Hobson (1931, pp. 444-53), we denote the Legendre function of complex 
degree P,,_;(s), usually referred to as a conal function, by 


K,(s) = P,-4(s). (5.9) 
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Che notation 
— s = 5. 0 
Ky (s) ds” P,,-4(8) (n 1, ) ( 


will also be used. 
The condition (5.4a) leads to the equations 


F(a, 27 +&,) K2(s)da = [s —t,]- —[s—cos (E, — (5.11a) 
0 


F(a, 4) KP (s)da = [s—t,] —[s—cos (5.116) 
0 
where t; = cosé;,4 = 1, 2. 

From the integral definition of the conal function K,,(s) given in Hobson (1931, 
p. 451), and the subsequent discussion given there, it can be shown that the 
expansion 


cosh — 77) 
s—tht=/2 K,(s)da 5.12 
is valid for 0 < € < 27. Further, by differentiation and a permissible mterchange 
of order of this with the indicated integration, we obtain 


[% cosh a(& 
cos 
V" coshan 


—7) da. (5.13) 
Both members of this are now multiplied by sin £ and the result integrated from 

&=£, tog = 27—(£,—{,). We obtain the representation 
[2 K2(s) 
— cos $_ ip) a 

|, (a? + 1) cosh az 
x {cos cosh —7) — cos — £,) cosh — £5 + §) 

—a{sin £, sinh —7)+sin (€,—§,) sinh (5.14) 


for 0 < £, < &. By an entirely similar procedure we find that 


oe) (1) 
[s—t,]-$—[s — cos (£2 —& = 4/2 


Jo (22+ 1) cosh am 
x {asin &, sinh a(&,—7) — asin (§,—&,) sinh a(&, —&)—7) 
— cos &, cosh —7) + cos cosh (5.15) 


for £, < & < 27. We may now replace the right-hand members of (5.11a) and 
(5.116), respectively, by (5.14) and (5.15), and we are then led to the equations 


F(a, £,+27) = sin (£, —,) sinh a(7 — £4 + 


+asin &, sinh «(&, —7) + cos —&,) cosh a(7 — + 
—cos, (5.16) 
for 0 < £, < &, and 
F(a, &) = {asin sinh «(£,—7) —asin sinh a(&,—&)—7) 


—cos £, cosh —7) + cos 5) cosh (5.17) 
for £, < < 2a. 
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From the second boundary condition (5.46) we obtain, assuming that (5.11a) 
and (5.116) are satisfied, 


K2(s {[s—t,]-* —[s—cos (5.18) 


1 
oF 


If (5.14) and (5.15) are now inserted in (5.18)—(5.19), the differentiation and 
integration which occurs in each of the resulting equations may be interchanged, 
and we obtain finally 


oF /2 
(a, + 27) {sin (5 — £,) cosh — + + sin cosh 20) 
for 0 < £, < &, and 
oF /2 
(a, = —*—— {sin £, cosh a(£, —7) —sin (£,—£) (5.21) 


cosh az 


If we let 
F(a, £) = cos coshaé, &) = cos€ sinh (5.22) 
Fy,(a,£) = sin§ coshaé, = sing sinha, 
and (cos cosh a(§ — 7) — cos — £5) cosh — +77) ) 
—a{sin sinh a(§ —7) —sin (€ — £5) sinh a(§ —£,+7)] 
(0<£<&),| 


ats 1) a(é 7) — cos (g &,) cosh = 7) 
—a[sin § sinh —7) —sin — £9) sinh a(£ — £5 —7)] 


(fo < § < 2m), 
(5.23) 


then (5.16) becomes, noting (5.7). 
F,,(a, + 277) + Fyo(a, + 277) B(x) + Foy (a, + 277) C(x) 
+ Fyo(%, §; + 277) D(x) = g(a, §)). (5.24) 
Further, (5.17), (5.20) and (5.21) yield 3 additional equations linear in A,..., D. 
These with (5.24) constitute four equations which yield unique solutions for 
A,...,D by the use of Cramer’s rule if the determinant 
| Fyy(a, + 27) + 27) 


| , , ‘ 


(5.25) 


is non-vanishing. After a rather extensive calculation we find that 


F (a, = a* sin? (£,—§1)— sinh? (Eo — 27) a, (5.26) 
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and it is easy to show that this is < 0ifa > 0. Fora = 0, F does vanish, however, 
and Cramer’s rule cannot be used, but we shall see later that this is a limiting 
case which causes no difficulty. The calculation of A, ..., D is nothing more than 
a formidable exercise in algebraic manipulation, and will be omitted. We give 
merely the result of inserting these quantities in (5.7): 


F cosh 
F(a,£) 


sin [A, sin (£,—£&) cosh a(£, — 27) — A, sin (£, — &) cosh — £)] 


a{(A, cosh a(7 — 7) — A, cos 7) sin (£, — €) sinh — &) 
+ (A, cosh — 277) — A, cos sin (£5 — sinh a (€; —§ + 277)] 
—sinh «(7 — 27) [cos —£)- £) sinh a(£, —§ + 27)]} 
— £)sinh a(£, —§ + 27) —T, sin (£, — §) sinh a(, 


+ sinh «(7 — 27) sin (£, sinh 


—asin7[I', sin (&, 


—T,sin (€,—&) sinh a(€, —§ + 27)], (5.27) 
where 7 = £,—£,, and 
A; = T;= (x €;) cosh an (¢= 1,2). 
(5.28) 


Here and elsewhere in this section the symbol ’ indicates differentiation with 
respect to 
As noted earlier, A, ..., D, and hence F, are not defined for « = 0. Examination 


of (5.27) reveals, however, that lim F(z, &) is finite, so that the singularity of the 
a—>0 
integrand of (5.6) at « = 0 is removable. 


6. Particular lenticular configurations of interest 
(a) The hemisphere 
In this case £, = $7, £5 = 37, & = 7, and we have 


(sinh? 37a — a) cosh am F(x, £) 


= eos ||? cosh — Q(x )sinh ac cosh 
a 


+1 
x sinh a(7—£)— N(x) lari cosh — £) + sinh 37a sinh — 
+sin& [sinh 37 sinh — €) + a? cosh «(32 — 


Nia | 
( cosh 37a — aQ(x) — cosh 377m sinh ina sinh a (377 — al, (6.1) | 

where M(x) = cosh 37a + «(sinh sinh $77), 
Nia 


)= 
= cosh ?7x + /2 cosh 


cosh }7a+asinh $72, (6.2) 
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It is hardly necessary to point out that much of the complexity (and little of the 
symmetry) of the general case is still present in this example. This is not sur- 
prising in view of the complex body shape considered. 


(b) The symmetrical (bi-convex) lens 


In this instance, if § = £, is one face of the lens, then the other is § = £ = 27 —§&,, 
and £, = 7. We find that (5.28) and (5.23) with these values give 


A, = A, = —2cosh cosh — 37) —asin&, sinh a(&, — 47)], 
= = 2cosh Jamsin £, cosh «(£, — 37), (6.3) 
F (a, = a? sin? 22, — sinh? 2aé,, (6.4) 


and, after some manipulation, the specialization of (5.27) for this case can be 
written in the form 


= a,£) = cosé cosh — 
+ £,) +a? sin €, sin cosh sinh ws, 


+sin sinh «(27 —&) [cosh £,(sin sinh + Da cos £,) 


[sinh a£,(—cos sinh 2a, 


+ a’ cos sin 22, sinh aé,]—T',‘¥ cos £, cosh ab], (6.5) 
where £,) = cos 2, — cosh (6.6) 
€,) = asin 2¢, —sinh 2a). 


Still further simplification can be achieved by noting that '¥ is a factor of F and of 
the coefficients of the A, terms. Removal of these common factors yields 


cosh am(x sin 2¢, + sinh F(a, &) 


= cos cosh a(27 — (asin &, cosh + cos, sinhaé,)—T, sin €, sinhag,| 


+sin & sinh «(27 — &) (a cos £, sinh «£, —sin cosh cos cosh 


(6.7) 


A further specialization of interest is obtained by setting £, = 37, which gives 
the sphere. In this case A, = A, = 0, and T, = —T, = 2cosh (az). It is easily 


checked that 
h —&) 
—9./94 208 
F(a, £) = 24/28 (6.8) 
and the integral in (5.6) then becomes 
) cosh — £) 
If we make use of the representation 
K (s) = coth an | (6.10) 
7 » (cosh u—s)# 
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and the discussion of pp. 451-3 of Hobson (1931), it can be shown that 


K(s) cosh — £) dx 


Differentiation of this and a permissible interchange with the integration on the 
right then yields 


2/2 0 cosh a7 
Thus (5.6) reduces to 
(=) 
st+t 
and, since (p/b)? = (s +¢)/(s—t), we obtain 
=- > (1 (6.13) 


which is the solution for the sphere given by Lamb (1945, p. 598). 


(c) The calotte, or spherical cap 


If £, = £,(= &,), the two bounding surfaces of the lens coincide, and the body 
becomes a portion of a spherical surface bounded by a circle of latitude. This is 
called a calotte, or spherical cap. As & goes from 47 to 7 the cap goes from a 
hemispherical surface of radius } to a flat disk of radius 6. Cases 0 < £) < 47are of 
less interest. An approximation to the former type of body is found in the cup of 
the meteorologist’s anemometer. 

In this case the procedure used in obtaining the equation (5.27) for F(a, &) is 
invalid. Nevertheless, this formula, with A, ..., Das determined there holds in the 
limit as £, > To show the expression [s —t]~? — [s — cos — is repre- 
sented by (5.15) for £, < § < 27 and by (5.14) with £, replaced by ¢,—2z for 
2n < £ < 2n+é,. In either case the resulting integral is combined with that 
appearing in (5.6). The result is a form of yy — is valid near the peak of the cap, 
i.e. in the neighbourhood of the point ( = £5, 7 = 0) at which the x-axis pierces 
the cap. 

We give the details of this process only for the case of the hemispherical cap. 
This is defined by £, = £, = £ = 47, and we have 

\2 
g(a, = gla, = {cosh am — x sinh jaz], (6.14) 
» 


= g'(a, £5) = cosh haz. 6.15 
(%,51) = 9 So) ( ) 


With the values of A; and I’; which these give, (5.27) then yields 


F(a,é) = {a2 cos & cosh a(7—& 
(%,€) (a2 + 1) cosh? am 
+ a[eosh am cos sinh a(37 —£)+sinh sin cosh «(37 — £)] 


+ cosh «($7 — £) [cosh cos — cosh sin £]}. (6.16) 
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We may write (5.6) for £, = 37 in the form 


J0 


F(a, &) (6.17) 


and making use of (5.6) and (5.14), respectively, as indicated in the preceding 
paragraph, 


(s —t)-4—(s—sin 
= in sinh —£) + cosh «(37 —&)] 
+ cos [a sinh (4a < < 27), 


(s —t)-* — (s —sin £)- 


=—,/2 {sin £[a sinh — £) — cosh «($7 — £)] 
+ cos sinh (27 < < 3m). (6.18) 

The insertion of these and (6.16) in (6.17) yields 

= — fa? cos cosh — + sinh «(42 —&) 

(a? +1) cosh az ‘ 


x [sinh }amcos§—acosh samsin£]}dx (6.19a) 
for $7 < ' < 27, and 


= 4Ur?(s—t)3 la 


— sinh — [sinh cos€+acosh Zamsiné]}da (6.196) 


s) 

cos cosh — 
1) cosham* 
for 27 < € < $7. 

Because of the formal nature of the operations yielding (6.19) it must be 
verified that the yy thus defined is actually the solution of the cap problem. This is 
found to be the case. It is clear that yy does not exhibit a singularity at the point 
where the x-axis pierces the cap. 


7. The flow about a pair of separated spheres 


The curves 7 = const. defined by the transformation (5.1) constitute a family 
of circles in the av-plane whose centres lie on the x-axis. All circles 7 > 0 lie 
entirely within x > 0 and enclose (b, 0); those for which 7 < 0 lie entirely within 
x < 0 and enclose (—6,0). Rotation of these curves about the x-axis yields a 
family of spheres similarly characterized by 9 <= 0 

We consider two spheres defined by 7 = 9, > 0 and y = 7, < 0. It is clear that 
neither lies in the interior of the other and that they have no points in common; 
we refer to them as separated spheres (figure 4). The region exterior to both is 
9, (0< § < 27). (7.1) 

In order to solve the Stokes flow problem for the region (7.1) we seek solutions 
of (2.9) which are periodic in & of period 27. We choose a decomposition for 7, in 
(3.4) to be a sum of (a), (b) and (c) of (3.2), so that 
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It is then easy to see that a suitable y is 
[A,e"" cosh 7 + B,,e-™" cosh 
n 
+C,e™sinh + D,,e-"" sinh 9] P?(cosé)} (7.3) 
which, by a redefinition of constants, can be written as 


\ n=1 


+ P® (cos tea 


Ficure 4. Separated spheres. 


In order to satisfy the boundary conditions on 7 = 7, and 7, we recall the well- 
known expansion formulas (Hobson 1931, p. 335) 


Jj2 P (t) > 0), (7.54) 

(s—t)4 
2 (7) < 0). (7.56) 

n=0 


We now differentiate (7.5a) with respect to ¢, multiply both sides of the result by 
sinh 7, and integrate with respect to 7 from 7, (> 0) to «©. In this way we obtain 


(s, —t)-* = =< > 9) (7.64) 


and in a similar way find that 


e(n— Dr) 


(s,—t)-# = p> P(t) < 9). (7.60) 


With the aid of these results we now see by reference to (7.4) that the condition 
that yy should vanish on 9 = 7, and 7, will be satisfied if we choose A,, ..., D,, to 
satisfy the relations 


n+3 n—4 
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The second boundary condition, that dy/on = 0 on » = 9, and 7, can be easily 
shown to reduce to 


n=Ni n=Ni 


n= 
(7.8) 
and hence will be satisfied if 
(n+1)A, (n—1) B, (n+ 1)C, — (n—-1) D, 
=(-1)'}/2 ini! — iil] (i = 1,2). (7.9) 


Equations (7.7) and (7.9) constitute 4 equations for the determination of the 
constants A,,,..., D,,,and hence the solution of the problem. We leave the problem 
at this point, since the solution has been given, although in somewhat different 
form, by Stimson & Jeffery (1926). 


8. The flow about a spheroid 


The solution of the Stokes flow problem for the prolate and oblate spheroids are 
implicit in the results of Oberbeck (1876) for the general ellipsoid. Since these 
problems furnish, however, striking examples of the ease with which problems can 
be solved by judicious choice of the representation formula for y, and since the 
results have a simpler form than that given by Oberbeck, the authors feel 
justified in sketching briefly the solution for these bodies by the methods of this 


paper. 


(b) 
Figure 5. (a) Prolate spheroid. (b) Oblate spheroid. 


The transformation (8.1) 


where z = x+ir and € = £+7y, serves to introduce elliptic co-ordinates in the 
ar-plane, line segments £ = £ = const., 0 < 9 < 7, of the plane being mapped 
into the upper halves (r > 0) of confocal ellipses in the 2r-plane, with foci at 
(+c¢,0). Rotation of such a curve about the a-axis generates a prolate spheroid 
(figure 5a) whose exterior is defined by 


<7). (8.2) 


We represent the solution of the differential equation L* ,y = 0 in the form 


= = yp). (8.3) 
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It is easy to see (Hobson 1931, p. 413) that Q,(s)P,(t) and Qi 8) Prt), 
n = 0,1,2,..., where s = cosh, t = cosy, are, respectively, functions 1 and 3 
which are aneliie in the region (8.2). Thus, at least formally, we may choose 


From the condition that dy/07 vanish on € = &, it follows at once that we must 
ane A,=0, A,,B,=0 (n> 2), 
and hence that 
= — Ay Qo(s) — By PP(t) Q2(s)}. (8.5) 


If we insert in this the expressions for P,(t) and Q,,(s), n = 0,1, and relabel the 
coefficients, the result is 


The boundary conditions (3.3) now serve to determine A and B, and we find that 


= 8(8§ — 1)/(s% + 1) —3(s2 + 1) In(s +.1)/(s— 1)| (8.7) 
(8+ 
where = cosh £, 
An easy calculation shows that the drag is 
-1 
P = 8np = 87 (8.8) 
The transformation 
z=csinh€ (c> 0), (8.9) 


yields oblate spheroids. The segments £ = £ = const., 0 < 9 < 7, of the ¢-plane 
are mapped into the halves of confocal ellipses with foci at (0, +c) lying in r > 0. 
The rotation of such a curve about the z-axis generates an oblate spheroid 
(figure 5b) whose axis coincides with that of 2. The exterior of the spheroid is 
(0< <n). (8.10) 
Again we represent the solution of (2.9) in the form 

= >). (8.11) 


In this case we find (Hobson 1931, p. 422) that P,(t)Q,(i7) and P®(t) QM(i7), 
n = 0,1,2,..., where 7 = sinh and ¢ = cosy, are, respectively, functions y and 
v3 which are regular in (8.10). We may thus assume y to have the form 


= 3Ur - Pal) + Bn (ir) (8.12) 
By a procedure similar to that used in the preceding case we find that the y which 
also satisfies the boundary conditions, and thus constitutes the solution of the 


problem is 7™(1 + 79)/(1 + 7) + (1 — 75) cot? 7) (8.13) 
t 2 \ To +(1 — 7?) To 


where 7, = sinh ). 
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The drag is found from (4.2) to be 
P = + (1 — 72) cot (8.14) 


The case of the flat circular disk is of particular interest, and is the special case 
of (8.12) obtained by setting 7, = 9. Accordingly, 


9 
and the drag is P = 16uU¢e. (8.16) 


It is not immediately obvious that (8.7) and (8.13) yield the same results as 
were obtained by Oberbeck, since he eschews the use of both the stream function 
and velocity potential, and deals directly with the velocity components, ex- 
pressing them in terms of the gravitational potential of the ellipsoid and related 
integrals. For the ellipsoid of revolution, these integrals reduce to types which 
can be easily evaluated and the authors have done this to verify that in such case 
Oberbeck’s formulas reduce to (8.8) and (8.14). The result (8.16) for the disk was 
also obtained by Oberbeck. 

We note, finally, that for c fixed and &, large, both spheroids closely approxi- 
mate a sphere of radius 4c e5o, and it can be shown, as one would expect, that both 
(8.8) and (8.14) become approximately 67U(4ce%), the drag of a sphere of 
radius $c 


9. The drag problem for the lens 


The drag of the axially symmetric body is given by (4.2). In the case of the 
lens-shaped body y, consists of the last two terms of (5.6), where F(a, ) is 
obtained from (5.27)—(5.28). The evaluation of the formidable expression which 
results we will not attempt. The bi-convex lens, because of symmetry, gives an 
integral which is considerably simpler, but still requires several numerical 
integrations which the authors have not carried out. There is, however, one 
special case of interest in which the drag can be computed,} viz. that of the 
spherical cap. From (6.17) we find that 


=4Ur +(s—t) j F(a, &) da (9.1) 
and using (5.1) and (2.2) we have 
where 
/2 
F(a, 27) ‘x? cosh am — a cosh am sinh tam + cosh? sam} (9.3) 
and K®(1) = —4(a?+}4) (9.4) 


+ Another such special case is the sphere, treated by Stokes. See §6 above, and the 
reference given there. 
35-2 
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(Neumann 1881). Thus we find that 


sinh 
cosham 


P = 4,/2mubU a*sech anda — 


2 Jo + coshamn 2] 

(9.5) 


These integrals are all available in Bierens de Haan (1939), except for the last one. 
It can be evaluated in closed form by a contour integration and subsequent 
summation of a series which we will omit. We then obtain for the drag 


da— sech am da 


1 2 
P= = 17: 525U by. 


This is somewhat in excess of the value 16Uby which was obtained for the flat 
disk in the last section. 


Body Drag 

Hemispherical cup 4 /2mubU (5 17-525ubU 
Flat disk 16nbU 
Sphere 674bU 


1 -1 
Prolate spheroid E +1) log 


Oblate spheroid 79) ¢ cot! + 
Two identical n(n+1) 4sinh? (n+ —(2n+1)*sinh?a 


sink 1- 
spheres 2sinh (2n + 1)a+ (2n+ 1)sinha 


TABLE 1 


Table 1 has been prepared to permit easy comparison of the cup drag with 
that of the other bodies which we have discussed. In every case 6 is the 
radius of the frontal area circle; sy and 7, are defined in §8, and a = cosh~!d, 
where d is the ratio of the distance between centres of the circles to their diameter. 
The result for the two spheres is taken from Stimson & Jeffery (1926). It wiil be 
noted that the drag of the hemispherical cup lies about midway between that of 
the flat disk and the sphere. Calculation also shows that the drag of the oblate 
spheroid always lies between that of the disk and sphere, and the drag of the 
prolate spheroid is always greater than that of the sphere. 


This work was done at the National Bureau of Standards — Dr Payne 
was engaged as a consultant there, and was supported by the the U.S. Air Force, 
through the Office of Scientific Research of the Air Research and Dev elopment 
Command. 
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t 
The influence of hole dimensions on static t 
pressure measurements 
By R. SHAW 
Department of Mechanical Engineering, University of Liverpool : 
(Received 1 July 1959) i 
The pressure measured at a static pressure hole differs slightly from the true d 
static pressure, by an amount which depends on the hole size and shape. The ” 
present investigation extends the range of previous work to determine the error 
in static pressure measurement in incompressible turbulent flow. The observed ” 
static pressure was always greater than the true static pressure. The results are 8) 
presented in dimensionless form as a function of the Reynolds number based on "7 
hole diameter and friction velocity. be 
w 
1. Introduction ne 
Previous investigations have shown that the observed static pressure is di 
influenced by the dimensions of the static pressure hole. The flow near to the wall pl 
is disturbed by the hole, and the fluid in the hole is set in motion by the fluid 
passing it. An infinitely small hole has been assumed to give the correct reading, wl 
and the experimental results for relative errors have been extrapolated to zero Ve 
hole diameter in order to assess the absolute error. m 
Investigations have been carried out by Fuhrmann (1912), Allen & Hooper TI 
(1932), Myadzu (1936), and Ray (1956) for incompressible flow and by Rayle di 
(1949) for both incompressible and compressible flow. Ray correlated his results mi 
by employing the principle of dimensional similarity and expressed the error as tes 
7 a function of the Reynolds number and the length/diameter ratio of the hole. | inc 
Rayle correlated his results on a Mach number basis. The present experiments Re 
extend the range of results for incompressible turbulent flow. 
Fuhrmann employed a model in a wind tunnel in which the air velocity was 
32 ft./sec. The model consisted of a body of revolution with very thin walls and 
three interior compartments. A 0-031 in. diameter hole located in the wall of | wh 
the centre compartment served as a reference, while two test holes located of 
: diametrically opposite at the front compartment were varied in size from C : 
: 0-005 in. diameter to 0-043 in. diameter. For sharp-edged holes, negative errors co" 
: were observed. These errors increased rapidly with increase of hole size up to ba 
: 0-016in. diameter, the error then becoming nearly constant. rec 
Allen & Hooper experimented with a 12in. diameter pipe which had been bel 
cleaned by scraping. The static pressure holes were located in }in. diameter ’ 
removable plugs which were scraped flush with the bore of the pipe. The experi- hol 
ments were carried out with water flowing at velocities from 4 to 7-3 ft./sec. For soli 


holes varying in diameter from j; to 33 in. no error was recorded. An jin. radius val 


‘ 
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on an 4 in. diameter hole gave a positive error, but smaller radii and 4, in.) gave 
the same reading as a sharp-edged hole and were therefore recommended in order 
to ensure freedom from burrs. 

Myadzu employed a channel of square section with water flowing at velocities 
up to 14 ft./sec. Adjustable plugs with static pressure holes varying from 0-004 to 
0-157 in. diameter were used to obtain the error relative to a line of reference holes 
equally spaced along the channel length to record the fall in pressure. The 
absolute error was found by graphical extrapolation, and the results indicated a 
positive error which increased linearly with hole size but was independent of 
velocity. The depth of hole was significant; the error was constant for length/ 
diameter ratios greater than 2, but decreased for smaller ratios, until below a 
ratio of about 0-4 negative errors were indicated. 

Rayle located a 1 in. diameter test section at various distances downstream of a 
nozzle, and used insert plugs with static pressure holes of length/diameter ratio 
greater than 2-5 and diameter varying from 0-006 to 0-125in. Both air at mean 
velocities from 400 to 900 ft./see and water at velocities from 22 to 31 ft./see were 
employed in order to observe the effects of compressibility. Positive errors were 
recorded, and these errors were found to increase with increase in hole diameter, 
with increase of Mach number and with reduction of the distance between the 
nozzle and test section. Various edge forms were tested for 0-032 and 0-046 in. 
diameter holes. A radius was found to increase the positive error, whilst a chamfer 
produced a negative change. 

Ray’s experiments were carried out using a sugar solution, the concentration of 
which was varied to give kinematic viscosities between 0-01 and 0-03 ft.?/sec. 
Velocities up to 12ft./sec were employed. The test section was rectangular, 
measuring 2-8 x 1-6in., with static pressure holes located in the longer lower side. 
The diameter of the holes was varied from 0-039 to 0-393in., and the length/ 
diameter ratio from 1-75 down to 0-1. The diameter of the connexion to the 
manometer was also varied, being smaller than the hole diameter for most of the 
tests. The results indicated that the error was positive and increased with 
increase of the Reynolds number and reduction of the length/diameter ratio. 
Ray’s results expressed in similar form to that used in the present report give 


where AP is the absolute static pressure error, 7, the wall shear stress, / the length 
of the hole of diameter d, v the kinematic viscosity, and p the density. Here 
C = f{l/d} varies from 0-5 when //d = 1-75 to 1-08 when l/d = 0-1. His results 
cover the range 1-7 < R < 31-6, where R = (d/v) ,/(79/p) is the Reynolds number 
based on hole diameter and friction velocity. Negative pressure errors were 
recorded only when the length/diameter ratio was very small and the enlargement 
behind the hole considerable. 

Thom & Apelt (1957) have calculated the pressure in a two-dimensional static 
hole at low Reynolds numbers, by means of an arithmetical method for the 
solution of the Navier-Stokes equations. The solution is for laminar flow and is 
valid for the Reynolds numbers u,w/v up to 1, where w, is the velocity at the 
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centre line and w the width of a two-dimensional static hole. It is suggested that 
there might be two different functions for the dimensionless pressure error, one 
for laminar flow and one for turbulent flow, connected by a range of transition 
values, as in the case of resistance to flow in a pipe. The solution is modified to 
give the dimensionless pressure error for the three-dimensional case of a circular 
hole. If Nikuradse’s results for the resistance to flow in pipes are extrapolated to 
very small Reynolds numbers, this solution can be expressed in similar form to 
that used in the present report. At the upper limit (u,w/v = 1) for which the 
solution is valid, AP/7, = 0-088 and (d/v) (7 9/p) = 1-4. 


2. True static pressure 


The flow behaviour at the static pressure hole is shown in figure 1. The deflexion 
of the streamline into the hole is confirmed by the analysis for a two-dimensional 
hole by Thom & Apelt (1957). The curvature of the streamline is such as to 
increase the static pressure in the hole above the true value for the pipe. There is 
also an eddy or system of eddies set up in the hole and a Pitot effect at the down- 
stream edge of the hole. These three factors combine to give the net pressure 
error. For a fixed speed of flow it is probable that these effects will decrease as the 
hole size decreases, so that AP > 0 asd —> 0. 


Ficure 1. Flow behaviour at the static pressure hole. 


3. Local dynamical similarity considerations 

We assume initially that the static hole is small compared with the diameter of 
the pipe (or thickness of the boundary layer on a plate), that the geometry of the 
hole and manometer connexion remains similar, and that the pipe surface is 
smooth. Under these circumstances the pressure error for a finite hole size can be 
influenced only by the hole diameter d, the fluid density p, the absolute viscosity /. 
and the local flow conditions at the surface. Now for a distance y up to approxi- 
mately one-fifth of the pipe radius or boundary-layer thickness from the wall, the 
local velocity wu is given by 


from which it is apparent that w at distance y from the wall depends on 7», p and jv. 
It is therefore the velocity gradient, described by the wall shear stress 7), which 
defines the local flow conditions. 
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Since for a very small hole the fluid is not disturbed by the hole, it may be 
assumed that as d tends to zero the measured pressure tends to the true static 
pressure. (See next section where this is discussed further.) Accounting for the 
variables suggested above, the static pressure error for a finite hole size can 
therefore be written 

AP = f{T.p. 4,4}, 


AP d 


Allowing for changing hole and connexion geometry and increase of hole 
diameter until the hole-diameter/pipe-diameter ratio is significant, equation (1) 
can be rewritten AP 1 D D,\ 

p’ ad’ ad’ 
where D is the diameter of the connexion to the manometer and D,, the diameter 
of pipe. 


(2) 


4. Behaviour of AP/7, as (d/v) (7)/p) > 9 

It is impossible to experiment with extremely small holes or with very small 
velocities, so that (d/v) ./(7)/p) cannot be reduced to indefinitely smail values. It 
is therefore impossible to establish the true static pressure directly, and extra- 
polation is necessary. This procedure could lead to inaccuracies unless there is 
some theoretical guide to the behaviour of the error as (d/v) ,/(T)/p) > 9. 

For any given fluid (density p and absolute viscosity j fixed) the friction 
velocity ./(79/p) falls as the flow velocity is reduced, and the region in which the 
viscous stresses are important increases in thickness. Close to the surface the 
viscous stresses are large compared with the dynamic stresses, since the velocity 
u— 0. The effect of density on the motion is therefore negligible, and it can 
therefore be eliminated from the variables involved. Thus AP = f{79, u,d}. Only 
one dimensionless group can be formed, so that 


AP/t,. = const. (3) 


Since it is inferred that AP/7, > 0 as d > 0, equation (1) shows that AP/7, must 
also approach zero as (d/v) ,/(7)/p) > 0. Thus the constant in equation (3) must 
be zero, and therefore AP/7, = 0 as (d/v) ./(79/p) becomes small. 


5. Behaviour of AP/7, when (d/v) ./(79/p) is large 

For a pipe of large diameter (or very thick boundary layer on a plate), in- 
creasing the static hole diameter for a fixed velocity will result in a greater flow 
into and out of the hole and the turbulent velocity component v’ will no longer 
vanish at the hole position. Thus, for large holes, turbulent stresses may outweigh 
the viscousstresses. Also the Pitot effect at the downstream edge of the hole, which 
is a dynamic effect, will outweigh the viscous effects. Hence the viscosity may 
now be eliminated from the variables. Thus AP = f {7 , p,d}, and AP/7, = const. 
as only one dimensionless group can be formed. Thus the curve of AP/T,) vs 
(d/v) ./(7/p) must be as shown in figure 2. 
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Dimensionless pressure 
error (AP/7,) 


Reynolds number (d/v) ./(7)/p) 


Figure 2. Predicted curve of dimensionless pressure error against Reynolds number. 


6. Apparatus 


The general assembly of the apparatus is shown in figure 3. Air was drawn 
through 30 ft. of 2in. internal diameter brass pipe by either a three-stage centri- 
fugal blower or a small single-stage blower, and the air flow was measured by a 
}-radius Pitot tube flow meter (Preston 1950). The pipe inlet was shaped to 
reduce losses and was followed by a transition ring to promote turbulence 
(Preston 1958). 

2 fi. in 


Pitot-static flow meter}, | | 


Butterfly valve in 


tin. I.D. blower 
delivery pipe 


Flexible tube 


Test section 


Static taps 


Pipe entry and 
transition ring 


Three-stage blower 


Induction motor 


FiaurE 3. General assembly of apparatus for measurement of static pressure. Internal 
diameter of pipe = 2in.; mean flow velocity 38-212 ft./sec.; all pipe joints carefully 
blended and flanges dowelled for accurate location. 


The test section shown in figure 4 was located at 130 diameters from the inlet. 
in a region of fully developed turbulent flow. It consisted of a 10in. long 2in. 
internal diameter pipe made from 2? in. diameter solid brass, split on a diameter 
to facilitate inspection of the static pressure holes. The test section and the 
immediate upstream pipe section were honed after assembly. All other pipe joints 
were made in such a way as to ensure no wall discontinuities occurred along the 
length of the pipe. 

The test section was provided with sets of static pressure holes at various axial 
locations. All holes were opened out in the boss to a connexion diameter equal to 
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twice the hole diameter, the length of hole being varied between tests to give the 
various length/diameter ratios. A Chattock gauge was employed to measure the 
small pressure differences, and the connexions to this gauge were made by rubber 
tubing of diameter equal to twice the hole diameter and length equal to 240 times 
the hole diameter. 


A 
Section AA 
ection 
6 f 
Allen 
screws] | 
gin. dia. || |" 
taper pins) 
, 
dia. hole B. dia. hole 
10 in. 
Static hole diameter (d) in inches 
Position 1 2 3 4 5 6 7 8 
AA 0:0635 0-0940 0:1890 0-0275 0-0635 0-1610 0:1275 0-0165 
BB 0-0635 0-0310 0:1560 — 0:1275 0-0635 0-0940 0-0635 


Static hole counterbore dia. D = 2d. 
Static hole length J varied. 


Ficure 4. Test section. etail of static holes 


7. Procedure 

An attempt was made initially to use insert plugs located in a boss brazed on to 
the 6ft. pipe section near the downstream flange, but it proved extremely 
difficult to hone the plug flush with the pipe wall, and the errors due to a pro- 
truding or recessed plug were of the same order as the errors due to hole size. The 
special test section was then machined and, in the first series of tests, was honed 
after drilling the holes. The test section was then split, and the appropriate drill 
shank inserted into each hole from the inside wall position to remove the internal 
burrs. The holes were then inspected for cleanliness and size using a traversing 
microscope. Honing and visual inspection were repeated until all the holes 
appeared to be free from burrs. These preliminary results were reported (Shaw 
1958), but subsequent tests suggested that the sharp-edged form of the hole was 
destroyed by repeated honing. A Talysurf surface measuring instrument was 
then used to check the wall profile immediately before and after each hole, and 
this showed that the repeated honing had slightly rounded the edge of the holes, 
the effect being most marked with the larger holes. Since Rayle had shown that 
a radius tends to increase the positive error, it was apparent that new tests were 
necessary and that even greater care would be required in preparing the holes. 

A new set of holes was located at section AA (figure 4), two reference holes 
0-0635 in. diameter being provided. In this instance the test section was carefully 
honed before the holes were drilled, until the surface roughness indicated by the 
Talysurf instrument was less than 50 micro-inches. The holes were then drilled and 


i 
: 
RY 
5 
wh 
Ti- 
ya 
to 
1ce 
d 
nal 
ly 
et. 
in. 
ter 

he 
its 
he 
ial 
t 


556 R. Shaw 


the burrs removed by careful hand polishing with fine grade emery paper wrapped 
round a semi-cylindrical wooden block, the surface finish and edge form being 
repeatedly checked on the Talysurf instrument until all burrs were removed. 
Internal burrs were removed using the drill shank and brass polish. Typical 


00015 
= 1 Section BB. I/d=1. 0-0310in. dia. hole (upstream edge) 
uv 
a 4 
4 
A 
0 1 4 1 J 
0-001 0-010 0.020 0-030 0-040 0-050 0-060 
¢ 4 Axial distance (in.) 
- ] Section BB. i/d=1. 0:0635in. dia. hole (upstream edge 
4 
| 
2 Te 
3 4 \ 
3 4 \ 
2 | 
0-001 0-010 0-020 0.030 0-040 0-050 0-060 
: Axial distance (in.) 
= Section BB. i/d=1. 0-0940in. dia. hole (upstream edge 
v 
3 
0001+ 0-010 0-020 0-030 0-040 0-050 0-060 
Axial distance (in.) 
a Section BB. I/d=1. 0:1275in. dia. hole (upstream edge 
3 
4 
0-001- 0-010 0-020 0-03¢ 0-040 0-050 0-060 
4 Axial distance (in.) 
] Section BB. i/d=1. 0:1560in. dia. hole (upstream edge 
g 4 
a L L L 1 j 
0-010 0-020 0-030 0-040 0-050 0-060 


Axial distance (in. 


Ficure 5. Talysurf surface measuring instrument records. 


Talysurf records are shown in figure 5; the apparent inclination of the holes is due 
to the form of the Talysurf stylus. The length of the holes was reduced between 
tests using end mills of appropriate diameter to give //d ratios of 6, 4, 2, 1 and 0-5, 
care being taken with the smaller holes to ensure that the inside wall of the pipe 
was not deformed. 

The experiments were repeated for length/diameter ratios of 2, 1, and 0-5 using 
a further set of holes located at section BB (figure 4), the previous holes having 
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been plugged and the test section honed. In this case two additional precautions 
were taken. First, the dimensions of the pipe bore were checked using a Mercer 
cylinder gauge to ensure freedom from taper and ovality. These checks indicated 
that the maximum taper was 0-0003in. per in. at section BB and that the 
maximum ovality was 0-0013in. Second, in order to facilitate the drilling and 
end milling of the various holes, semi-cylindrical blocks were made to reinforce 
the pipe wall during these operations. 

In the initial tests a null-reading inclined-tube type of micromanometer was 
used, but the results suggested that different viscous damping effects were pro- 
duced by the asymmetry of the limbs. The Chattock gauge was therefore used 
because of its symmetry, but the dimensions of the manometer connexions were 
taken proportional to the hole size to preserve dimensional similarity up to the 
manometer. A further inconsistency in the results persisted at low velocities where 
the errors appeared to be comparatively large. Since this was thought to be due 
to some spurious effect originating at the three-stage blower, which was severely 
throttled at low flows, some of the later tests were carried out using the small 
single-stage blower. This eliminated the discrepancies. 

For each test hole, readings of pressure error, relative to a reference hole, were 
recorded against the dynamic pressure at the ?-radius position. The duplicate 
reference holes were used to assess the pressure difference between two nominally 
similar holes, and in general these differences were found to be small. At the 
maximum velocity the average difference was 0-007 in. of water (i.e. the absolute 
pressure error for the 0-0635in. diameter reference hole, when J/d > 1-5, is 
0-0360 + 0-0035 in. of water). The maximum difference between reference holes 
was 0-013in. of water. Oscillations of the Chattock gauge bubble made accurate 
readings difficult at some velocities, the maximum bubble fluctuation being 
equivalent to about +0-005in. of water. All observed relative errors were 
corrected to errors relative to the mean of the two reference holes. 

On several occasions during the tests, readings were taken of pressure at the 
test section and also of the pressure drop over a 4ft. length of test section. The 
mean velocity of flow was calculated from the observed dynamic pressure at the 
3-radius position, assuming the flow to be incompressible. The flow meter was 
calibrated initially by traversing one of the four Pitot tubes. The wall shear stress 
was calculated from the observed pressure drop at the test section, account being 
taken of the fluid acceleration due to compressibility, and the resistance coef- 
ficient y = 7,/}pu?, was plotted logarithmically against the pipe Reynolds number 
D,,u,,|v, where u,, is the mean flow velocity. These results were found to agree 
with Nikuradse’s results for turbulent flow in a smooth pipe. 

The results for length/diameter ratios between 1-5 and 6 were found to be the 
same, within the limits of experimental accuracy, and were therefore plotted on 
the one graph. After cross-plotting, graphs were drawn showing the pressure 
error, relative to the 0-0635 in. diameter reference holes, against diameter of hole 
for various mean velocities and length/diameter ratios (e.g. figure 6), and these 
curves were extrapolated to zero hole diameter to obtain absolute pressure 
errors. The extrapolations of these curves to zero hole diameter were made to 
fulfil two requirements, namely (1) that the curves pass as closely as possible 
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through the points obtained for small diameter holes; and (2) that the intercepts 
(at d = 0) are such that the absolute pressure errors, thereby obtained, give 
results which are consistent with the requirement, from similarity considerations, 
that the dimensionless pressure error is a function only of the Reynolds number 
(d/v) ,/(79/p). Graphs of dimensionless pressure error against Reynolds number for 
various length/diameter ratios are shown in figures 7, 8 and 9. 

0060 
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0-100 
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x 
J 


— 0-020 


—0-040 


FiaureE 6. Pressure error (relative to 0-0635 in. diameter reference holes) vs diameter of 
hole for various mean velocities: l/d ratios 1-5 to 6. 
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FIGURE 7. Dimensionless pressure error against Reynolds number for l/d ratios 1-5 to 6. 
Hole diameter (in.): @, 0-175; +, 0-150; A, 0-125; m, 0-100; 0, 0-075; 2, 0-0635; x , 0-050; 
A, 0-025. 
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Finally, four 0-0635in. diameter holes with length/diameter ratios of 4 were 
drilled at a new axial location to assess the pressure error due to burrs projecting 
into the pipe. To obtained a variety of burrs, the holes were drilled with the semi- 
cylindrical support blocks in place, with the same drill speed, but with different 


A, 20 
be 
ry 
is 

2 
= 
WE 
tt 
3) 
= 

x 
0 100 200 300 400 500 600 700 800 


Reynolds number (d/v) 4/(79/p) 


FicurE 8. Dimensionless pressure error against Reynolds number for //d ratio of 1. Hole 
diameter (in.): @, 0-175; +, 0-150; A, 0-125; m, 0-100; 0, 0-075; x , 0-050; A, 0-025. 
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FicuRE 9. Dimensionless pressure error against Reynolds number for //d ratio of 0-5. 
Hole diameter (in.): @, 0°175; +, 0-150; A, 0-125; g, 0-100. 


feed rates; larger burrs being produced by higher feed rates. The burrs produced 
were measured on the Talysurf instrument and found to be 0-0003, 0-0004, 
00017 and 0-0022in. in height. The hole with the 0-0003in. burr was then 
polished by hand until the burr was removed. A set of readings was then obtained 


= 
ts 
he 
e 
Ss, Ab 
+ 
he 
asi 
‘ 
6. 
st? 


560 R. Shaw 
for the pressure error of the holes with drill burrs relative to this square edge hole f 
for various mean velocities. The burrs were then reduced by hand polishing to t 
0-0001, 0-0004 and 0-0009 in. height, respectively, and the procedure repeated. 
00015 g 
= Reference hole (downstream edge li 
3 
2 r 
if 
0 4 J 
0.001 - 0-010 0.020 0-030 0-040 0-050 
Axial distance (in.) 
= ] Hole with 0-0001 in. burr (downstream edge) 
s | 
4 
0001" 0-010 0-020 0-030 0-040 0-050 
> Axial distance (in.) 
c 
] Hole with 0:0004 in. burr (downstream edge) 
a 
3 4 
3 
4 
0-010 0-020 0-030 0-040 0-050 
Axial distance (in.) 
Ss | Hole with 0-0009 in. burr 
(downstream edge) 
a 
a 
0 A. i 4. i 
0-002" 0-010 0-020 0-030 0-040 0-050 
> 7 Axial distance (in.) 
= Hole with 0:0017 in. burr (downstream edge) 
3s 
a 4 
4 
0-010 0-020 0-030 0-040 0-050 
Axial distance (in.) 

FicureE 10. Talysurf surface measuring instrument records. the 
lalysurf records showing the shape of the burrs are shown in figure 10, and the ae 

. . . . 

results are plotted in dimensionless form in figure 11, together with the hole size D 
error curve which is replotted on this curve for comparison purposes. be 

8. Discussion of results and comparison with results of other workers dri 
The dimensionless curves of AP/7, against (d/v) ./(7)/p) (figures 7-9) show the | ghc 
form which was predicted by the dimensional analysis and the physical con- Sin 
siderations. For //d > 1-5, the error is independent of the //d ratio (figure 7), but poi 
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for //d = 1-0 (figure 8) and //d = 0-5 (figure 9) there is a progressive reduction in 
the error. 

Inspection of figure 7 shows that the points for d = 0-175in. lie below the 
general curve. This might suggest that the hole diameter d is then sufficiently 
large compared with the pipe diameter D,, for the dimensionless group D,,/d to be 
significant. However, figures 8 and 9, for length/diameter ratios of 1 and 0-5, 
respectively, do not exhibit the same behaviour. A study of figure 6 indicates that 


if the measured relative error for the 0-1275 in. diameter holes is slightly excessive, 
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Dimensionless pressure error (AP,/T,) or (AP/T) 


Reynolds number (d/v) ,/(79/p) 


FicureE 11. Dimensionless pressure error against Reynolds number for 0-0635 in. diameter 
holes with varying heights of radial burr: ratio 4. —, AP./Ty); AP/T). 


then the curves for each mean velocity would not tend to flatten at the larger hole 
diameters as much as is shown, and the points on figure 7 would then lie more 
nearly on a unique curve. It would therefore appear that the dimensionless group 
D,/d is not significant for values greater than 10. Experiments with very large 
holes would be necessary to establish the onset of a D,,/d effect. 

Figure 11 shows the effect of burrs formed in a natural manner during the 
drilling of the hole. For comparison, the hole-size error of a well-finished hole is 
shown, and it is immediately apparent that the effect of burrs is very significant. 
Similarly, specks of dust may have an important effect. From a practical stand- 
point, the smallest hole may not necessarily have the least error, since for a given 
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height of burr or speck of dust the effect will be larger for the small hole than for 
the larger one, and may well outweigh the error due to hole size. 

Although no evidence of negative errors was obtained in the present experi- 
ments it is believed that such errors could be produced by holes of small length/ 
diameter ratio which open out to very large diameters. 
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Dimensionless pressure error (AP/T,) 
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Figure 12. Dimensionless pressure error against Reynolds number for J/d ratio 1-5. 
Comparison with Ray’s results (extrapolated for R > 31). (a) Ray’s results for 1/d = 1:5, 
connexion diameter < d. (b) Present experimental results for //d = 1-5, connexion 
diameter = 2d. 


Figure 12, showing Ray’s results (extrapolated for R > 31) indicates that the 
dimensionless pressure errors in his case are considerably greater than the values 
found in this series of experiments. Some difference in the results can perhaps 
be attributed to the different connexion diameters employed, but the main 
source of difference occurs because of the method of extrapolation to zero hole 
diameter. In the first instance, Ray did not experiment with holes smaller than 
0-039 in. diameter, and secondly, the intercepts for d = 0 were not selected in the 
way described above. Ray’s correlation was carried out by assuming that the 
dimensionless error was proportional to R~” (where n > 0) which was approxi- 
mately true for his experimental points, but only for d > 0-039 in. 

Figure 13 gives a comparison with Rayle’s results for incompressible flow. 
Rayle plotted a dimensionless pressure error AP/4pu?, against hole size for Mach 
numbers M of 0, 0-4 and 0-8. In his calculations for incompressible flow, the 
dimensionless errors for velocities varying from 22 to 31 ft./sec were averaged. 
Reynolds number effects being neglected. This gave a unique curve for M = 0. 
The variation of error with velocity was not clearly shown, because the velocity 
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range was limited, and an inclined manometer was employed to measure relative 
pressure errors. 

Since it is the wall shear stress 7, which influences the error, it is desirable to 
plot the dimensionless pressure error as AP/7,. Therefore, taking Rayle’s results 
for water flow in the case when the holes are located 27 diameters downstream of 
the nozzle and assuming the flow to be fully developed (Nikuradse suggests that 
for turbulent flow the ‘length of transition’ is 25-40 diameters), and reworking 
Rayle’s results, using y = 7,/4pu?, from Nikuradse’s curves for smooth pipes, we 
obtain the graph shown in figure 13. All the points are for relatively high mean 
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FicurE 13. Dimensionless pressure error against Reynolds number for J/d ratio > 1-5. 
Comparison with Rayle’s results for incompressible flow. 


velocities of flow, but the hole size varies from 0-006 to 0-125in. diameter. If the 
flow in Rayle’s experiments is not fully developed, then the correct 


= H[du/dy], 


will be greater than assumed, and this should reduce Rayle’s AP/7, to give even 
closer agreement than is indicated in figure 13. 


9. Conclusions 

For static pressure holes which are normal to the pipe wall, square-edged, and 
of length/diameter ratio between 0-5 and 6, with manometer connexion 
diameters twice the hole diameter, the pressure error is positive in turbulent flow 
and is a function of the Reynolds number (d/v)/,/(79/p) only, as shown in figures 7 
to 9, for Reynolds numbers up to 800, The dimensionless error tends to zero for 
small Reynolds numbers, but increases progressively more rapidly with Reynolds 
number up to about 300, and then progressively less rapidly up to Reynolds 
number of 800, when the rate of increase is quite small. The error is only influenced 
by the length/diameter ratio of the hole for values of this ratio below 1-5, and is 
then found to be smaller. 
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For a smooth pipe with incompressible turbulent flow and a static pressure hole 
diameter one-tenth of the pipe diameter, the static pressure error reaches about 
1°, of the mean dynamic pressure at a pipe Reynolds number of 2 x 10°. 

With a ,}, in. diameter hole of length/diameter ratio 4, the error due to a drill 
burr, which projects into the main stream, is approximately equal to the error due 
to hole size for a burr height ¢ of 0-0005in. (e/d = 1/127), and is approximately 
seven times the error due to hole size for a burr height of 0-0020 in. (e/d = 1/31-7), 


see figure 11. 


The author wishes to thank Prof. J. H. Preston for permission to use the 
facilities of the Department of Fluid Mechanics and to acknowledge the assistance 
given by him during frequent discussions on this project. 
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’), (Received 29 July 1959) 

- This article has a twofold purpose: (1) to analyse the available theoretical and 

we experimental knowledge concerning flow in the inlet region of a smooth round 
tube, and (2) to point out that the e® amplification factor method apparently 
predicts natural transition correctly over a significant fraction of the entire inlet 
length of the tube. The successful prediction indicates, but does not prove, that 

54 flow in a smooth round tube becomes turbulent at higher Reynolds numbers 


because transition occurs in the inlet length—not in the fully developed Poiseuille 
régime. The close agreement between theory and a test result obtained by Pfen- 
ninger indicates that the e® method is valid for a wide variety of flows having 
tic x Reynolds numbers of transition ranging from 570,000 to 40 million. The results 
are applicable to both plane and axially symmetric flows. 


1. Introduction 

ent In 1956 Smith & Gamberoni showed by analysis of a large number of experi- 
ne. ments that boundary-layer transition will occur when Tollmien-Schlichting 
aaa waves reach an apparent amplification ratio of about e®. One experiment, not 


treated in the original work, was that performed by Pfenninger (1950, 1951a, b) 
at the Northrop Aircraft Co., in which he carefully measured the transition point 
ow ina long (59 ft.) tube of 2 in. diameter containing a flow with very low turbulence. 
This type of test is of particular interest because the problem of transition in 
a round tube has been paradoxical for many years. Theory predicts that fully 
developed parabolic flow is stable with respect to small rotationally symmetric 
disturbances at all Reynolds numbers; yet the flow is known to become turbulent 
at rather low values of the Reynolds number. The reasons Pfenninger’s tube-test 
was not included in the original study were two: (1) the flow was not of a boundary- 
layer type except very near the entrance, and (2) no amplification rate chart 
applicable to this kind of flow was in existence. Subsequently, efforts have been 
made to analyse the flow, and several observations concerning the state of 
knowledge about this problem were made that are believed to be of general 
interest. This paper then has a twofold purpose: to report that the e® factor does 
seem to hold for the tube flow studied and, secondly, to call attention to and 
discuss the unsatisfactory state of knowledge concerning this classical yet still 
paradoxical flow problem. 
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2. Experimental data 


_Many tests have been made to learn the tube Reynolds number Ry = Udlv 
(U = mean velocity, d = diameter, v = kinematic viscosity) at which the flow t 
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For a smooth pipe with incompressible turbulent flow and a static pressure hole 
diameter one-tenth of the pipe diameter, the static pressure error reaches about 
1% of the mean dynamic pressure at a pipe Reynolds number of 2 x 10°. 

With a ;/; in. diameter hole of length/diameter ratio 4, the error due to a drill 
burr, which projects into the main stream, is approximately equal to the error due 
to hole size for a burr height ¢ of 0-0005in. (e/d = 1/127), and is approximately 
seven times the error due to hole size for a burr height of 0-0020 in. (e/d = 1/31-7), 
see figure 11. 


The author wishes to thank Prof. J. H. Preston for permission to use the 
facilities of the Department of Fluid Mechanics and to acknowledge the assistance 
given by him during frequent discussions on this project. 


REFERENCES 


ALLEN, C. M. & Hooper, L. J. 1932 Piezometer investigations. Trans. A.S.M.E. 54, 
HYD 54-1. 

FUHRMANN, G. 1912 Theoretische und Experimentelle Untersuchungen an Ballon- 
modellen. Dissertation, Géttingen. 

Myapzu, A. 1936 Influence of type and dimensions of pressure hole on the recorded static 
pressure. Ingen.-Arch. 7, 35. 

Preston, J. H. 1950 The } radius Pitot tube flow meter. Engineer, 27 October, p. 400. 

Preston, J. H. 1958 The minimum Reynolds number for a turbulent boundary layer and 
the selection of a transition device. J. Fluid Mech. 3, 373. 

Ray, A. K. 1956 On the effect of orifice size on static pressure reading at different 
Reynolds numbers. Ingen.-Arch. 24 (3), 171. Translation in A.R.C. Rep. no. T.P. 498. 

Raye, R. E. Jr. 1949 An investigation of the influence of orifice geometry on static 
pressure measurements. M.S. Thesis, M.I.T. 

SHaw, R. 1958 The measurement of static pressure. Report Advisory Group for Aero- 
nautical Research and Development, no. 163, p. 1. 

Tuom, A. & Apert, C. J. 1957 The pressure in a two-dimensional static hole at low 
Reynolds numbers. A.R.C. Rep. R. & M. no. 3090. 


T 

e) 

ti 

fl 

be 

ré 

ni 

ar 

1. 

Wi 

tre 

at 

in 

Tl 

dis 

at 

We 

] lay 

ap 

mé 

kn 

int 

See 

al dis 

pa 


Ow 


Remarks on transition in a round tube 
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This article has a twofold purpose: (1) to analyse the available theoretical and 
experimental knowledge concerning flow in the inlet region of a smooth round 
tube, and (2) to point out that the e® amplification factor method apparently 
predicts natural transition correctly over a significant fraction of the entire inlet 
length of the tube. The successful prediction indicates, but does not prove, that 
flow in a smooth round tube becomes turbulent at higher Reynolds numbers 
because transition occurs in the inlet length—not in the fully developed Poiseuille 
régime. The close agreement between theory and a test result obtained by Pfen- 
ninger indicates that the e? method is valid for a wide variety of flows having 
x Reynolds numbers of transition ranging from 570,000 to 40 million. The results 
are applicable to both plane and axially symmetric flows. 
1. Introduction 

In 1956 Smith & Gamberoni showed by analysis of a large number of experi- 
ments that boundary-layer transition will occur when Tollmien-Schlichting 
waves reach an apparent amplification ratio of about e®. One experiment, not 
treated in the original work, was that performed by Pfenninger (1950, 1951 a, b) 
at the Northrop Aircraft Co., in which he carefully measured the transition point 
in a long (59 ft.) tube of 2 in. diameter containing a flow with very low turbulence. 
This type of test is of particular interest because the problem of transition in 
a round tube has been paradoxical for many years. Theory predicts that fully 
developed parabolic flow is stable with respect to small rotationally symmetric 
disturbances at all Reynolds numbers; yet the flow is known to become turbulent 
at rather low values of the Reynolds number. The reasons Pfenninger’s tube-test 
was not included in the original study were two: (1) the flow was not of a boundary- 
layer type except very near the entrance, and (2) no amplification rate chart 
applicable to this kind of flow was in existence. Subsequently, efforts have been 
made to analyse the flow, and several observations concerning the state of 
knowledge about this problem were made that are believed to be of general 
interest. This paper then has a twofold purpose: to report that the e® factor does 
seem to hold for the tube flow studied and, secondly, to call attention to and 
discuss the unsatisfactory state of knowledge concerning this classical yet still 
paradoxical flow problem. 


2. Experimental data 


_Many tests have been made to learn the tube Reynolds number R; = Udlv 
(U = mean velocity, d = diameter, v = kinematic viscosity) at which the flow 
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becomes turbulent; for example, see Goldstein (1938), Schlichting (1955), or 
Prandtl & Tietjens (1934). Yet very little study has been made of the nature of 
the flow in the initial portion of the tube, the inlet length, and one wishing to 
study the problem in detail will find that a complete set of reliable test data does 
not exist. Three of the principal sets of data will be mentioned here. One is due 
to Nikuradse (see Prandtl & Tietjens 1934). It is complete in a static sense, 
consisting of a full set of velocity profiles covering the entire inlet length. 
However, all the results are presented in two small graphs with neither explana- 
tion nor description of the test apparatus. Therefore the data have a low degree of 
accuracy, if for no other reason than the small size of the figures. 

Another set is Pfenninger’s data. In this case great care was used in taking 
measurements, the apparatus was very good, and the scale was large. Neverthe- 
less, the data are not entirely satisfactory, because the tests were in the nature of 
spot checks rather than a complete, systematic survey. For example, transition 
location is reported only for one tube Reynolds number. No curve of R,, (the 
Reynolds number based on 2) for transition versus R, is given. Several velocity 
profiles were measured, but these were taken under different conditions and in 
slightly different apparatus from those for which the transition measurement was 
made. In addition, the Northrop data cover only a small fraction of the entire 
inlet length. 

Recently a third set of data have become available: the results of tests run by 
Reshotko (1958) in an extremely fine experimental arrangement located under- 
ground in a nearly isothermal environment. But again the data cannot be called 
complete. The purpose of the tests was to study the stability of fully developed 
Poiseuille flow. Consequently, only the downstream portions of the inlet region 
were measured, and the Reynolds number could not be increased sufficiently to 
induce transition in the inlet length. (Another experimental investigation has 
been performed by Leite (1959), its purpose and interest being quite similar to 
Reshotko’s. However, his data will not be cited here, because Reshotko’s appara- 
tus was considerably superior, and because Leite covered only the fully developed 
flow régime. Nevertheless, his results supply a great deal of information on the 
stability of fully developed flow.) 

Figure 1 compares transverse velocity profiles according to Nikuradse, 
Pfenninger and Reshotko. In order to prepare this figure, it was first necessary to 
cross-plot Pfenninger’s and Reshotko’s results after the fashion of figure 13 of 
Prandtl & Tietjens (1934). The profiles so found are seen to have none too good 
agreement, differing by as much as 9 % at the centre. Notice also that Pfennin- 
ger’s data cover only the first part of the inlet length, whereas Reshotko’s cover 
only the last part; and since there is no overlap, they cannot be compared 
directly. The entire length corresponds to a value of x/aR,, of about 0-26. (Because 
x/aR,, has been the most common length parameter in the past, it is used here: F, 
is the tube Reynolds number Ua/v, where a is the radius—not the diameter—and 
U is the mean velocity.) The question of whether the differences are attributable 
to graph-reading accuracy, experimental accuracy, differences in entrance condi- 
tions, or other factors cannot be answered. Some of the same data are shown in 
a different and perhaps clearer fashion in figure 2. In figure 3 the core velocities 
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U, are compared. Again, Nikuradse’s and Pfenninger’s data are seen to be in 
only fair agreement. Reshotko’s data are mostly in good agreement with 
Nikuradse’s, but they too fail to overlap Pfenninger’s. This figure shows clearly 
that Pfenninger’s tests covered only a small fraction of the entire inlet length. 

In summary, only three sets of detailed measurements of mean velocities exist. 
One, Nikuradse’s, is complete but of uncertain accuracy; the other two are of 
good accuracy but incomplete for the present interest. 
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FicurE 1. Measured and calculated laminar velocity profiles at several positions in the 
inlet portion of a round tube. Experimental: —-_J—O—(, Nikuradse; —O—O—O, 
Pfenninger (195la, Fig. 3); —A—A—A, Reshotko. Theoretical: —-—-——, Lang- 
haar; — - —, Punnis; - - - - - , Tatsumi. 


3. Recent theoretical studies 


Between theory (Lin 1955; Corcos & Sellars 1959) and experiment (Leite 1959), 
it has become fairly well established that the asymptotic (parabolic) flow in 
a round tube is stable to rotationally symmetric small disturbances at all Rey- 
nolds numbers. Consequently, if stability theory is to explain natural transition 
for flow of low turbulence in a smooth round tube, it must explain it either by 
consideration of the inlet region alone or else by consideration of disturbances 
that lack rotational symmetry, in which case the Poiseuille régime may also enter 
the picture. Since it is known that transition can occur in the inlet length, recent 
efforts at explanation have dealt with flow in this portion. 

In order to apply stability theory, accurate primary velocity profiles must be 
available. The other results of the older theoretical attempts (Goldstein 1938) to 
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Fiacure 2. Comparison of measured velocity profiles with calculations by Falkner-Skan 


piecewise method. —O—O—O, Pfenninger (195la, Fig. 3); —D—O—O, Nikuradse ; 
- \— A— A, Reshotko. 
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FiGuRE 3. Comparison of several measured and calculated values of velocity in the core of 
a round tube. —[(J—, Nikuradse; —O—, Pfenninger (19516, Fig. 5, U, 9 = 26-28 m/s). 
A, R,= 2050, A R,=3800, A R,=18,000, Reshotko. --- - - , Tatsumi; — -—, Punnis; 
— — —, Langhaar. 
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find them are unsuitable, and we shall consider only the three most recent 
attempts to solve the problem. Those to be considered are by Langhaar (1942), 
Tatsumi (1952), and Punnis (1947). Results of these attempts are presented in 
figures 1 and 3. It is seen that the theoretical values are not only in disagreement 
among themselves but also in disagreement with the two sets of experimental 
results. For many purposes the agreement indicated might be considered good, 
but for purposes of stability calculations the accuracy is inadequate. 

The only calculations of the stability of these inlet profiles are those by 
Tatsumi, who computed the neutral stability loops for several values of x/aR,. 
His computed values of R,+,, (= Ud*/v, where 6* is the displacement thickness 
and U is the edge velocity) are shown in figure 4. In Pfenninger’s test the measured 
location of transition was at x/aR, = 0-017. Hence Tatsumi’s stability calcula- 
tions did not cover a sufficient portion of the inlet region to make possible cal- 
culations of wave amplification, even if he had computed complete stability 
maps—which he did not. 


4. A rough calculation of the neutral-stability and transition loci 


By now it should be obvious why Pfenninger’s tube-test was not included in the 
original correlation studies. Since both experimental and theoretical data 
needed for amplification studies of inlet flow were found wanting, a decision was 
made to apply ordinary methods of boundary-layer calculation to learn whether 
or not they might produce a useful result. The method chosen was exactly that 
used in the earlier studies, the Falkner-Skan piecewise method (Smith 19565), 
using Mangler’s transformation. The boundary-layer profiles resulting from the 
calculation are characterized by Hartree’s £. Once the variation of # along the 
axis of the tube is established, the stability or degree of instability of the flow can 
be established from Pretsch’s set of amplification rate curves (see Smith & 
Gamberoni 1956). 

The tube was treated as one of decreasing diameter, whose effective area 
decreased inversely as the core velocity supplied by the experiment. Figure 2 
shows the boundary-layer profiles calculated by the piecewise method. Both this 
figure and figure 1 show that the profiles have about the same level of accuracy at 
z/aR,, = 0-03 as those computed by the more involved methods of calculation. 

Tatsumi calculated stability of the flow only for values of w/aR, < 0-004. In 
this range the piecewise method has very high accuracy and is clearly much more 
rapid than Tatsumi’s method. 

It was not known whether the location of the neutral-stability curve was very 
sensitive to the shape of the core-velocity distribution. To learn the answer, 
neutral stability for two-dimensional disturbances was calculated using two dif- 
ferent core-velocity distributions, Pfenninger’s and Nikuradse’s. The results are 
shown in figure 4. The two different distributions cause appreciable difference in 
the calculated stability of the flow. 

How does the stability calculated indirectly from the core-velocity-distribu- 
tion data agree with that calculated directly from experimental velocity profiles? 
Pfenninger measured carefully the velocity profiles at a station 39-64 ft. from the 
entrance. Because he tested at several velocities, the profiles cover a small range 
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of values of x/aR,. The neutral stability of four of these profiles covering the 
x/aR, range were calculated by means of Lin’s short formula. The necessary 
velocity and derivatives were obtained by fitting a 20-point least-squares cubic 
to the boundary layer in the range 0 < y/d* < 1. The results, plotted in figure 4, 
show that the direct calculation of stability from the experimentally measured 
boundary-layer profiles produces considerably higher values of R,s,, than those 
determined indirectly from the core-velocity distribution. Furthermore, beyond 
x/aR,, of about 0-015 the stability of the experimental velocity profiles undergoes 
a rapid increase, while that determined from the core-velocity distribution begins 
to decrease. 
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Fiaure 4. Variation of neutral stability along the inlet length of a round tube according to 
several calculations. (a) Experimental, direct computation from measured velocity profiles 
(Pfenninger 195la, Fig. 3) using Lin’s short formula. The flow is treated as if it were two- 
dimensional. (b) Pfenninger data (1951, Fig. 5, U,,9 = 26-28 m/s) calculated by piecewise 
method (Smith 1956) and Pretsch charts (Smith & Gamberoni 1956). (c) Nikuradse data, 
calculated by piecewise method and Pretsch charts. (d) Tatsumi. 


Reshotko’s measurements are in general agreement with these other data but 
his velocity traverses included too few points for the difficult problem of com- 
puting stability of an experimentally measured velocity profile. 

Tatsumi’s calculation is also shown in figure 4. It represents the stability of 
rotationally symmetric disturbances, however, whereas the other data represent 
the stability of two-dimensional disturbances. For its interest, the critical 
Reynolds number of two-dimensional Poiseuille flow is also included. This value 
has been obtained from the commonly accepted value of R = 5300, where R is 
based on the half-width of the channel and the velocity at the centre of the channel. 
With respect to the other data this value is surprisingly low. 
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Figure 5 shows the point of transition calculated by the e® method, the 
neutral-stability locus according to Tatsumi, and the locus according to the 
present method of calculation. For purposes of calculating transition to the very 
end of the inlet length, Nikuradse’s core-velocity distribution was used. Also to 
be seen on the figure is the experimentally measured transition point supplied by 
Pfenninger. The calculated point nearby is computed in the same way as the full 
curve based on Nikuradse’s data, except that it -ses the Pfenninger core-velocity 
distribution existing during the test when the tiansition point was observed. The 
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Figure 5. Transition and neutral-stability loci according to several calculations. Also 
shown is Pfenninger’s transition measurement. (a) Calculated, using data from Pfenninger 
(19516) Uj) = 26-28 m/s. (6) Experimental Pfenninger (19516) U, , = 26-28 ms. (c) e® 
transition curve, Nikuradse core velocity. (d) Neutral curve, Nikuradse core velocity. 
(e) End of inlet length. (f) Neutral curve, Tatsumi. 


most interesting feature of this figure is the close agreement of the predicted and 
the experimental transition point. Moreover, the small difference between the 
location of the circled Pfenninger point and the Nikuradse locus indicates that the 
differences in core-velocity distribution do not cause significant changes in the 
calculated location of transition, even though the differences in neutral-stability 
locus were appreciable (figure 4). Since theory and Pfenninger’s test are in good 
agreement, and since boundary-layer calculations gain in accuracy near the 
entrance, it is probable that the curve of #, for transition versus z/aR, is fairly 
accurate for values of z/aR, less than about 0-02. At higher values of 2/aR,, 
the predicted transition curve should rapidly develop major errors. Clearly, 
this portion of the curve is very conservative, as can be realized by the fol- 
lowing consideration. Farther forward in the inlet, the flow is accelerating, and 
as a result it is a flow having positive values of Hartree’s / (or Pohlhausen’s A). 
Downstream of the inlet length, where the flow has become parabolic, the 
velocity gradient becomes zero. Obviously, boundary-layer calculations would 
predict an ordinary Blasius boundary layer of low stability for this region. But 
from the correct calculations we know that the flow in this region has infinite 
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stability, at least for axially symmetric disturbances. In spite of the ultra- 
conservative nature of the calculations for the larger values of z/aR,, the curve 
shows that at R,, = 7-1 x 10° transition will occur at the very end of the inlet 
length. 


5. Miscellaneous remarks on the transition phenomenon 

If all disturbances are damped in the parabolic region+ and if the flow succeeds 
in reaching this region in the laminar state, it may remain laminar forever. 
Conversely, if natural transition does occur in the tube, it will occur in the inlet 
portion. In view of the very conservative nature of the calculations for 
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Friaure 6. General nature of the transition locus in the 2/aR, vs R, plane. (a) Transition 
locus (this part can be calculated with acceptable accuracy). (b) More nearly correct 
result. (c) Parabolic flow. (d) Conservative boundary-layer calculation. 


x/aR,, > 0-02, it is clear from the figure that laminar flow can exist to the end of 
the inlet length, at least to R, = 7-1 x 10°, that is, to R, = 14-2 x 108. Experiment 
has shown that the flow can remain fully laminar to values of R, of about 30,000. 

The more nearly correct relations between FR, and x/aR, for transition should 
be qualitatively as sketched in figure 6. Beyond the inlet, 2, for transition 
becomes constant, and #, for transition becomes infinite. 

In the R,vs Ux/v plane the transition locus should look somewhat like the 
curve in figure 7. Although this figure is principally presented for schematic 
purposes, it is actually calculated from the upper line of figure 6, and in fact the 
curve to the right of 2, = 40,000 represents the same values as those shown in 
figure 5, that is, the portion to the right of Pfenninger’s point is the theoretical ¢? 
prediction. Pfenninger’s measurement has been added for orientation purposes. 
Clearly, if R, is very great, the flow begins as ordinary Blasius flow, because the 


+ This is a controversial assumption, because the flow has been proved stable only for 
rotationally symmetric disturbances. Little is known of the stability with respect to other 
types of disturbance. 


| 

be 
be 

TI 

to 
tur 
ap] 
tul 
str 
rea 
; the 
floy 
| nw 
tra 
tur 
bet 
enc 
: its 

ciel 
a fi 


t 


or 
er 


Transition in a round tube 573 


boundary-layer thickness is so negligible compared to the pipe radius that the 
core velocity is constant initially. Basically the transition locus varies in a hyper- 
bolic fashion and approaches a vertical asymptote, which is R, ..:+.- 

The above facts and considerations indicate the following sequence of events. 
The flow begins as a very thin boundary layer. At first it is so thin that it is stable 
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FicuRE 7. Nature of the transition locus in the R, vs FP, plane. 
@. Pfenninger’s measurement. 


to small disturbances. But shortly, as it grows, it becomes unstable and dis- 
turbances begin to amplify. If the flow is such that the tube Reynolds number is 
greater than the critical value, the disturbances grow until they undergo an 
apparent amplification ratio of about e®, at which point transition occurs. If the 
tube Reynolds number is reduced, the point of transition moves farther down- 
stream, either in z measure or 2/aR, measure. Finally, a Reynolds number will be 
reached at which an amplification ratio of magnitude e® is just attained. Below 
the value e®, waves will grow for a while, but then they will enter the more stable 
flow farther downstream and will die out. An experiment in which Reynolds 
number is varied by reducing the velocity in a particular tube would show the 
transition point moving downstream slowly at first. But when sufficient dis- 
turbance amplification could no longer be reached, the point of demarcation 
between the laminar and turbulent regions would wash downstream to the very 
end of the tube. The value of 2, for transition would then have a discontinuity in 
its variation with R,,. 

Experiments seem to indicate a discontinuity, but they have not been suffi- 
ciently thorough to settle the question. Whether the discontinuity could truly 
exist cannot be answered by theory because damping rates are known to be 
a function not only of disturbance amplitude but also of wave length. Distur- 
bances other than two-dimensional or rotationally symmetric ones would further 
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modify the picture, and the process just described is by no means the only one by 
which turbulent flow can develop. 

The old question of whether or not the inlet portions of the flow may be 
neglected can be made to appear absurd by a little consideration. Consider flow 
in a tube of 2in. diameter (the same as Pfenninger’s). The inlet length ends at 
z/aR,, of about 0-26, say 0-24 to make the arithmetic easier. Then for a 2 in. tube, 
the end of the inlet is at x = 0-24aR, orx = (0:24 + 12) R, = 0-02 R,. If R, = 104, 
x = 200ft. The first 200 ft. of any laminar flow can hardly be ignored! 


6. The e® correlation plot 

In Pfenninger’s test the turbulence was very low and the tube was very smooth. 
Thus his measurement qualifies for inclusion in the correlation plot of Smith & 
Gamberoni (1956). This plot compares measured transition points with ones 
calculated by assuming transition occurs when Tollmien-Schlichting waves are 
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Fricure 8. Transition data, correlation between measurement and prediction by stability 
theory. Data from Smith & Gamberoni (1956) except for Pfenninger’s measurement. 


Standard deviation: Type of body No. of points % deviation 
2-dimensional 31 11-5 
3-dimensional 10 17-2 


Note: (1) Solid points represent flight-test data. (2) Circled points represent bodies of 
revolution. (3) Open points represent wind-tunnel data. (F,)tr = (Utr2tr)/v. 
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amplified by a factor of e®. The result is shown in figure 8. The numerical values 
are (R,),, = 40-4x 10°, with U, = 35-17m/sec, = 18m, v = 15-65 x 
m?/sec. The Reynolds number (R,),,, based on momentum thickness, is approxi- 
mately 3600. The remaining points can be identified by reference to the earlier 
work. Considering the rough method of calculation, the agreement between 
measurement and prediction is, to say the least, a remarkable coincidence. Here 
R, is based on the local core velocity U,, not U. The additional point increases 
the range of validity of the correlation appreciably, from a transition value of 
R,, of about 17 million to over 40 million. The correlation now covers a range of 
Reynolds numbers #, for transition from about 570,000 to slightly over 40 million. 
Both extremes are axially symmetric flow, the lower being flow about a sphere, 
the upper being flow inside a round tube. Standard deviations for the sets of 
points are included. The values are not consistent with those shown in the earlier 
work, because there was an error in the original calculations. 


7. Conclusions 


1. Existing experimental data for the flow in a round tube are inadequate and 
unsatisfactory. 

2. Existing theoretical solutions for the velocity distribution in the inlet of 
around tube are of insufficient accuracy, at least for the requirements of stability 
theory. 

3. Provided the boundary-layer thickness is less than about half the tube 
radius, conventional boundary-layer calculations can predict velocity profiles 
along the inlet of a round tube with about the same accuracy as the more elaborate 
methods. 

4. The existence of turbulent flow in a smooth round tube can be explained by 
consideration of the inlet region. Transition can readily occur in this portion 
because the velocity profiles have a finite stability limit. 

5. Apparently the e® method will predict the location of transition in the inlet 
length of a smooth round tube throughout a significant fraction of its length. 

6. At transition, the apparent amplification factor for Tollmien-Schlichting 
waves is found to be a constant over a very large range of Reynolds numbers and 
variety of flows, from R, = 570,000 to over 40 million. The constant factor is 
about e°, regardless of body shape. 


REFERENCES 

Corcos, G. M. & Setuars, J. R. 1959 J. Fluid Mech. 5, 97. 

GoLpsTEIN, S. (Ed.) 1938 Modern Developments in Fluid Dynamics. Vol. 1, Ch. VII. 
Oxford University Press. 

LANGHAAR, H. L. 1942 Trans. A.S.M.E. 64, A-55. 

Leite, R. J. 1959 J. Fluid Mech. 5, 81. 

Lin, C. C. 1955 Theory of Hydrodynamic Stability. Cambridge University Press. 

PFENNINGER, W. 1950 Experiments with Laminar Flow in a Two-inch-diameter 40-foot-long 
Tube at high Reynolds Numbers. Northrop Aircraft Co. Rep. no. AM-128. 

PFENNINGER, W. 195la Further Laminar-flow Experiments in a 40-foot-long 2-inch- 
Diameter Tube. Northrop Aircraft Co. Rep. no. AM-133. 


be 
at 
4 
04, 
& 
are 
‘ a 


576 A. M. O. Smith 


PFENNINGER, W. 19516 Further Laminar flow Experiments in a Tube at high Reynolds 
Numbers. Northrop Aircraft Co. Rep. no, AM-147. 

Pranpti, L. & Trersens, O. G. 1934 Applied Hydro- and Aeromechanics. New York: 
McGraw-Hill. 

Punnis, B. 1947 Zur Berechnung der Laminaren Einlaufstromung im Rohr. Bericht 
47 P/03 Max Planck Inst. fiir StrOmungsforschung. 

ResHoTKo, E. 1958 Experimental Study of the Stability of Pipe Flow. 1. Establishment of 
an axially-symmetric Poiseuille Flow. Progress Rep. no. 20-364 Jet Propulsion 
Laboratory, Pasadena, California. 

Scuuicutine, H. 1955 Boundary Layer Theory. New York: McGraw-Hill. 

Surru, A. M. O. 1956a Transition, Pressure Gradient and Stability Theory. Proc. IX. Int. 
Congress of Appl. Mech., Brussels. 

Smitru, A. M. O. 1956) J. Aero. Sci. 23, 901. 

Smrru, A. M. O. & GamBERONI, N. 1956 Transition, Pressure Gradient and Stability Theory. 
Douglas Aircraft Co. Rep. no. ES 26388. 

Tatsumi, T. 1952 J. Phys. Soc. Japan, 7, 489. 


an 


| 
| 
| 
Tl 
4 ar 
sti 
ar 
co 
4 ve 
| us 
| so) 
th 
en 
th 
on 
no 
mé 
be 
(1s 
pe 
me 
diz 
pu 
Th 
an 
re\ 
gas 
fra 


The linearized flow of a dissociating gas 


By J. FK. CLARKE 
College of Aeronautics, Cranfield, Bucks 


(Received 1 August 1959) 


The equations for planar two-dimensional steady flow of an ideal dissociating gas 
are linearized, assuming small disturbances to a free stream in chemical 
equilibrium. 

As an example of their solution, the flow past a sharp corner in a supersonic 
stream is evaluated and the variations of flow properties in the relaxation zone 
are found. Numerical illustrations are provided using an ‘oxygen-like’ ideal gas 
and comparisons made with a characteristics solution. The flow past a sharp 
corner can be studied in a conventional shock tube and it may be possible to 
verify the present theory experimentally. In particular it may prove feasible to 
use the results to obtain a measure of the reaction rates in the gas mixture. 


1. Introduction 


When the chemical composition of a gas changes by chemical reaction a further 
source of dissipation is present in the flow field. The reactions are natural, 
thermodynamically irreversible, processes and as such lead to the production of 
entropy. The important difference between this type of entropy production and 
that associated with the transport phenomena is that it does not depend explicitly 
on the gradients of velocity, temperature and concentration. Thus its influence is 
not necessarily confined to the interior of boundary layers or shock waves, but 
may spread over the entire flow field. 

The equations governing the flow of a chemically reacting gas mixture have 
been derived previously, for example, Kirkwood & Wood (1957), and Clarke 
(1958a). In the present paper an attempt is made to consider flows with small 
perturbations, dealing only with the two-dimensional steady problem. The treat- 
ment is simplified by considering the dissociation reaction in a symmetrical 
diatomic gas A,, each A, molecule being made up from two A, atoms. In the 


pure gas the reaction is i 
A,+ A324, + Ay, (1) 


The species A, can be either A, or A,, either is assumed to be equally effective, 
and k, and k, are the overall specific reaction rate constants for the forward and 
reverse reactions. 

Further simplification is obtained by assuming that A, is an ideal dissociating 
gas (Lighthill 1957). In that case, writing c, for the equilibrium atom mass 
fraction of A, atoms in the mixture, we have 


¢2|(1—c2) = (pg RT |pW,) exp (— DW,/RT) (2) 
37 Fluid Mech. 7 
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(R is the universal gas constant, 7’ the absolute temperature, p the pressure, W, 
the molecular weight of the molecules and D the dissociation energy per unit 
mass). pq is a characteristic dissociation density which, for the ideal gas, is 
assumed to be constant. 

The flow equations are linearized, assuming small perturbations to the mean 
flow. Explicit solutions for the variation of pressure, density, atom concentration 
and temperature on the wall behind a sharp corner in supersonic flow are found as 
an example. It is found that the zone of influence of the corner is bounded 
upstream by the ‘frozen flow’ Mach line and also that the flow is not of the simple 
wave type as it would be in an ordinary inert gas flow. 


2. Two-dimensional steady flow 

In the particular case of two-dimensional steady flow, with velocity com- 
ponents wu and v directed along the x- and y-axes, respectively, the equations for 
a dissociating diatomic gas become* 


Cu ou Op 
(3) 
4 
| 
c)—c?} = 0, (5) 
+ pa {K(1 —c}=0, (6) 
os 
Tu 1—c)—c?} (4, = 0. (7) 


Equations (3) and (4) are the momentum equations (p is the mixture density) and 
equation (5) is the continuity equation for the atomic species (mass fraction c). 
T is a characteristic reaction time, defined in terms of the specific reaction rate 
constant for recombination, k,, as 


T = W3/4k,p?(1+c), (8) 
(k, is measured in (mole/unit volume) per unit time). The quantity K is defined as 
K = (Wyk,/4pk,), (9) 


and can be related to the atom mass fraction under equilibrium conditions in some 
circumstances. Forexample, if we choose toevaluate an composition 
at local values of p and 7’ then it can be shown that AK = (1 +c)c?/(1—c?), ¢ being 
the actual concentration of atoms under these conditions (Clarke 1958a). This 
form will be found useful below, but is discarded in the linear theory in favour of 
evaluation of c, at local pressure p and entropy s. 


* A brief account of the equations (5), (6) and (7) is given in Appendix A. 
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Equation (6) is a rearranged form of the overall mass conservation equation 
(Kirkwood & Wood 1957). o is a function of the thermodynamic variables only. 
In the present case 


pp, (2) 1 (2) 
, 10 
Cy, 7 p\ec pT (10) 
op oh 
where ; C= (sr) (10a 


(f, is the volume expansion coefficient and C,, the specific heat at constant 
pressure, both for the mixture in a chemically frozen state.) 
The thermal and caloric equations of state for the ideal dissociating gas are 


p = +e) (R/W,) 7, (11) 
h = (44+c)(R/W,) T (12) 


respectively. 4 is the specific enthalpy and D is the dissociation energy per unit 
mass.* o can be evaluated via equations (10a), (11), (12). In the present case 


o = [(DW,/RT)+ 1] (13) 
Equation (7) is the entropy equation, derived from the energy equation 


0, (14) 


and the thermodynamic equation 
T ds = dh—p™ dp fz) de. (15) 


#4, and 1, are the chemical potentials of atoms and molecules, respectively. It can 
be shown (Clarke 1958a) that 
RT, 
= W, log , (15a) 
if c, is evaluated at the local p and 7’ values. 


Two important results follow from the above set of equations. First, 
eliminating p between equations (3), (4) and (14) shows that 


h+4(u?+v?) = const. = hy (16) 
along streamlines, as in the inert gas-flow case. Secondly, defining the vorticity 
= ov ou 17 

Oa ey’ 
equations (3), (4) and (15) show that 
ch es 


* Other authors have made use of the dissociation energy per molecule, written as 
d usually. The relation between D above and d is D = N,d/W, where N, is Avogadro’s 
number. The quantities W,D/R = d/k (k = Boltzmann’s constant) have the dimensions 
of temperature and are sometimes written as 7',. For oxygen T', = 59,000 °K. 
37-2 
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The operator is equal to (u/q) — (v/q) 0/ex, where q? = u? + v?, and denotes 
differentiation normal to the streamlines. Equation (18) is Crocco’s theorem in 
two dimensions, generalized to include the case of chemical reactions within the 
flow field (see Hayes & Wu 1958). hy will be assumed constant everywhere (the 
flow is assumed to originate in a region of constant stagnation enthalpy). Using 
equation (15a), equation (18) in these circumstances becomes 
os RT 1—c?) dc 
(ut = + (=) (19) 
It follows that €is only zero when the flow is in complete chemical equilibrium 
(c = c,) or when the flow is chemically frozen (¢ = const.). (Note that equations (7) 
and (15a) show that s = const. everywhere in these cases.) 
Suppose now that the actual atom concentration c differs but little from the 
local equilibrium value evaluated at local p and 7’, i.e. put 


C=C,+C, C<€C,. 


If we assume in addition that ¢ < 1—c,, the entropy equation (equation (7)) with 


the aid of equation (15a) can be written approximately as 
os 0s 
0 


ux 
Cx 
where 7, = W3(1—c,)/8k,p?¢,. (p, is the density which would occur at the local 
values of p and 7 if the local atom concentration was equal to the equilibrium 
value c,.) This result indicates that the entropy rise along streamlines is of second 
order in the deviation of concentration from its equilibrium value at local values 
of p and 7’. We now make a slight but decisive change in the interpretation of the 
local equilibrium state. It is assumed that c, from now on refers to an equilibrium 
composition at the local pressure and entropy values, and that c differs from this 
c, by an amount c’, i.e. we put ees (21) 

By analogy with the case for which c, is evaluated at local p and 7’, we now 
write the chemical reaction rate term as 


(1/7) {K(1—c)—c?} ~ —c'/7’, (22) 


thereby defining 7’. (It does not seem possible to evaluate 7’ a priori in the same 
way as with 7, above, since with c, evaluated at local p and s, K in equation (9) no 
longer becomes a simple function of concentrations.) Then equations (5) and (6) 
can be written in approximate form, 


+0, +8 += = 9, (23) 
Ox Oy oy Tf 
Cu ev c’ 
paz —pazoS = 0, (24) 
Ox oy Oy, T 
and by comparison with equation (20) we write 
0s R 
Ox oy 
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Since c, now refers to equilibrium at local p and s values, 


dc, oc, 
deca (=) de (26) 
Writing in future (ec,/Op), = A, (26a) 
for brevity, equations (25) and (26) give 
dc, 0, Op op R c’? (Cc, 
} Equation (23) can now be written approximately as 


since the term in c’?/7’ is of a smaller order of magnitude than c’/7’. (We show 
below that (0c,/¢s),, is less than order unity.) 
In order to find the quantity A we proceed as follows. We can write 


a= +(e), (28) 


where the subscript e is added to the last partial derivative to emphasize that it is 
to be evaluated with c following its equilibrium composition. Under conditions 
of chemical equilibrium ~, = “7, and equation (15) becomes 


T ds = dh—p-dp. 
At constant entropy, therefore, 


oh ah 
=] dp— = 0 


A can now be found via equations (28), (29), (11), (12) and (2). For the ideal 
dissociating gas equation (2) shows that 


(i 2p oT 


where we now write D’ for W,D/RT’, and equation (12) shows that 


RT (=) 


ah’ R RT, (2c, 
(=r) = (4+¢,) (D 1) 


p 


The final result for the ideal gas is therefore 


{ | 
\(D’ + + 


pa = ¢,(1—¢2) (30) 
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In order to find (dc,/6s),, we first write it in the form 


(=) 


whence it follows, since (¢s/07'),, = (eh/eT),,./T from the thermodynamic 


equation above, that ~ - 
08}, oT}, 
From the results given above, the value of the derivative in the ideal dissociating 
gas 1s R +1) (32) 
W,\ 0s}, 2(4+¢,)+(D’ +1)? ¢,(1 


It can be seen that R(0c,/és),,/W, is roughly of order D’-?. For an oxygen-like 
ideal gas D’ = 59,000/7' so that D’ is roughly of order 10 in the interesting range 
of dissociation. 


3. Small disturbance approximations 


The assumption will now be made that disturbances to the free stream are small. 
The free stream is taken to be a gas whose velocity is U, parallel to the z-axis, and 
whose pressure, density and temperature are given by p,,, p,, and T.,. The degree 
of dissociation in the free stream isc,, and this is necessarily an equilibrium value. 

Furthermore, we assume that c, differs but little from c,,, an assumption which 
is justifiable provided both p and s do not change appreciably from their free- 
stream values. The (dimensionless) disturbance quantities are defined as follows, 


Pp = p.(1+p’'), 

T = +7"), 

p =p..(1+p’), 

u = U(1+u’), } (33) 
v= Uv’, 

¢=¢,+¢, 


where all the primed quantities are very much less than unity. Choosing some 
suitable characteristic length* L, and writing 


y = Ia, (34) 
the flow equations can now be expressed in dimensionless form. In writing down 


these equations we assume that squares and products of disturbance quantities 
are negligible to a first order of approximation. Equations (3) and (4) become 


Ou Op’ 

ye 5 
dv’ Cp’ 

36 

Pa Px 0 (36) 


* There is no geometric length characteristic of the flow round a sharp corner. L is 
entirely arbitrary and is introduced here solely for convenience so that the equations can 
be written in dimensionless form. 
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Since A will not differ markedly from its free-stream value, equation (27) becomes 


op’ oc’ 
Paro ag 0, (37) 


where suffix 00 indicates evaluation at free-stream conditions and we write 
= +, U/L, (38) 
where 7,, is taken as the mean value of 7’, i.e 
T. = (39) 


(I is the ratio of a characteristic reaction time to a characteristic flow time.) 
Likewise equation (24) reduces approximately to 


Op’ ou’ Ov’ 
Po + — Tap = 9. (40) 


We may remark here that a,, the frozen sound speed, is defined as 


(suffixes s and cindicating that entropy and concentration are held constant during 
the differentiation). It readily follows from equation (15), which may also be 


T ds = de +pd(p) (1, — de 
(e is the specific internal energy = 3(R/W,) T +cD for the ideal gas), that* 


4+ep 


Differentiating equation (40) with respect to 7, the terms in 0?p’/0£ dy and 
0?u’' /of Oy can be eliminated in terms of derivatives of v’ by using equations (35) 
and (36). The result is 


dy’ 


ec’ 
U2 2 + 08? + On? 


— on = (42) 
Thec’ term in equation (42) can be eliminated by using equation (37) differentiated 
with respect to 7 and the result is 
0 a? | 2 2 
* At constant s and c the thermodynamic equation gives 
[1—p(eh/ép)7,.] dp = 
dp = 
But for the ideal gas h is not a function of p and e is not a function of p under present 
circumstances. The result given in equation (41) follows at once. 
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But it has been shown by Clarke (1958a) that 


where a,,, is the equilibrium speed of sound in the free stream. Defining the frozen 
and equilibrium Mach numbers 
M, = M, = U/a,., (45) 
it follows that equation (43) can be rewritten as 
In order to find the equation satisfied by the pressure perturbation, the process of 
elimination of the variables can now be repeated in precisely the same way as that 
outlined above, only this time a start is made by differentiating equation (40) with 
respect to £, and it then readily follows that p’ satisfies an equation identical in all 
respects with equation (46). 
The two equations will now be used to find the flow behind a sharp corner in 
a supersonic stream. Before doing so, however, it is interesting to note that in the 
limiting cases of [ = 0 and [' = «, equation (46) reduces to the ordinary wave 
equation for sound propagation at the equilibrium and frozen sound speeds, 
respectively. In practice 0 < [ < , and it is with these conditions that we will 
deal below. 


+(1—M?) (46)* 


4. Supersonic flow round a sharp corner 


The flow is assumed to be supersonic in the sense that both M, and M, are 
greater than unity. Figure 1 is a sketch of the general configuration. At point O 
the flow turns through a small angle — 0, so that the equations for v’ and p’ derived 


Gas in chemical 
equilibrium 


Region influenced 
U by the corner 


Figure 1. Supersonic flow round a sharp corner. 


above can be assumed valid descriptions of the field. The flow from £ = —« to 
£ = 0 is uniform and quantities there will be denoted by suffix oo. 
The boundary condition at the wall is given by 


vy’ = —tan&(l+w’), when 7» = —étand, 
but to the present order of approximation this can be replaced by 
when 7~ 0. (47) 


* The writer’s attention has been drawn to a recent paper by Moore & Gibson (1959) in 
which a similar equation is derived by a different approach. These authors solve a similar 
problem to that treated below by approximating equation (46) in the form of the telegraph 
equation via a suitable co-ordinate transformation. 
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Since the flow is supersonic, all the disturbance quantities are identically zero 
upstream of the corner* and, accordingly, we define the Laplace transform of a 
disturbance quantity by 


Ble, 9) = [oe nyexp (—28) dé, (48) 


(where 6’ is any disturbance variable). 
It follows that 0 satisfies the equation 


i ~(2B,)? | 0, (49) 
where z is the Laplace operator (see equation (48)) and p satisfies a similar 
equation. In equation (49) 

Bj = Mj-1, B2= Mi-1, B* = > 1. (50) 
We choose the following solution for @, 


v = const. exp ( — 2B, Se 


since it represents an outgoing wave motion. The transformed boundary condition, 
equation (47), is 


v(z,0) = —O/z, 
2 
whence v = —(9/z)exp (-=2, 4 (51) 


Since equation (49) is also satisfied by p we can write (where A is a constant to 


p = Aexp 


The transformed version of equation (36), however, is 


and it readily follows from equation (51) that 


Defining a pressure coefficient C, as 
Cy = 2PaP' [Pa (53) 
the value of this quantity on the wall (y = 0 to a sufficient order of accuracy) can 
be found from A 20 | nats | (54) 
2B,| 


* See remarks made in §5 below. 
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586 J. F. Clarke 


The transform on the right-hand side of equation (54) can be inverted (Erdélyi, 
Magnus, Oberhettinger & Tricomi 1954) and we find that 


Cow = — (20/B,;) lexp[ 1) — 1) 


+ | —4(B2+1) Wyaw). (55) 
0 


I, is the zero-order modified Bessel function of the first kind. Figure 2 shows 
(B,C,,,,/20) +1 plotted against x/7,,U(= &/T) for a typical value of B?, namely 
1-5. With this value for B?, M, = 2-12 and M, = 1-83, respectively, if 

= 1:35. 
This value of the speeds of sound ratio can arise, for example, at a pressure of one 
atmosphere and temperature of 4,250°K in an ‘oxygen-like’ ideal dissociating 
gas (see Clarke 1958a). 


0-20 - 
“ 
+ 
Asymptote 
owl 
0 04 08 1-2 1-6 20 


FicureE 2. Variation of pressure on the wall; B? = 1-5. © = characteristics solution. 


It can be seen from figure 2 that, in contrast to the inert gas Prandtl—-Meyer 
flow round a corner, the pressure is not constant on the wall downstream of the 
point O, but rises steadily as x increases. It can be shown that C,,,, > — 20/B, as 
£/I'> oo. The integral in equation (55) can be rewritten* as 


Ba exp[—}(B?-+1)W] aw. 
g/T 


When é/T > « the asymptotic form of J, can be used both in the first term of 
equation (55) and in the integral above (i.e. [,() ~ (27x)-4exp a). The asymptotic 
form of C,,,, is then given by 


— B,C /20 ~ [n(B?- 1) + B41 (56) 


where erfc is the complementary error function. The result quoted above follows 
in the limit as £/[' > oo, and (B,C,,,,/20)+1 tends asymptotically to the value 
0-183 in figure 2. The percentage variation of pressure behind the corner is quite 
significant therefore. 

As remarked in §2, the flow is not irrotational for 0< I < oo. It seems 
reasonable to suppose that the continuous pressure rise behind the corner results 


* Note: exp dg = (a? — 
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from the reflexion of the primary (expansive) disturbances from the vortex sheets 
as compression waves. The vorticity arises as a direct consequence of the (thermo- 
dynamically irreversible) chemical reactions. 

In order to find the density variations in the flow it is observed that the mass 
conservation equation in its more familiar form, namely 


in the present approximation. Using equation (35) it follows that 

Asymptote 


FicureE 3. Variation of density on the wall; B? = 1-5, (a,.,/a,.~)% = 1:35. 
© = characteristics solution. 


The transform of the density increment, /, is therefore (using equation (58)) 
p = 
In particular it is found from equations (54) and (51) that, on the wall, 


The first term here is simply (31?) times the pressure coefficient transform. 
Noting that 1+Tz is the transform of I'-!exp(—£/I) and making use of the 
convolution theorem, it can eventually be shown (see Appendix B) that 


B,pi,|OM} = exp —(B*-+1) £/2P 1) 


The value of —(B,p/,/9M?)—1 is shown plotted against 2/7, U in figure 3 for 
B? = 1-5 and (a,.,/a,,,)” = 1-35. By similar methods to those used for C,,,, we can 


show that p),> —0M?/B, as z/7,,U >. It is observed that the variation of a 
density behind the corner is smaller than that of the pressure in the present case 
and that p,,, decreases continuously along the wall. 4 


The values of atom concentration on the wall can be found as follows. Equa- 
tion (37) can be rewritten in the form 


0 
(c’ exp £/T) = — exp (€/T) 
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and it follows after an integration by parts that c’ on the wall is given by 
0 


Making use of the expression for C,,,, (equation (55)) the expression for c,, reduces 
to (see Appendix B), 


Cy = (p.. U?A,,0/B,) exp [ — (B? + 1) £/21] H[(B? 1) £/21). (59) 
The asymptotic value of c;,, as £/[’ > oo is therefore 
Cy ~ (P. U?A,,6/B,) [7(B?— 1) g/l exp (—£/P), (60) 


showing, as we might expect, that c;,, > 0 and the concentration approaches a new 
equilibrium value. 
The ‘small disturbance’ version of equation (27a) can be written as 


where it can now be verified (using equation (59)) that the last term there can be 
neglected in a first approximation, at least near the wall. (There seems to be no 
reason why it should not be negligible everywhere, but equation (59) verifies that 
it is so only near the wall.) Then 


Coun = = tPo U7AQC, (61) 


pw 


From the definitions of c, and c’ (note c—c,, = c,+c’) and using the results 
equations (60) and (61), it now follows that 


= —(p,A,Mj O(4+¢,,)/3B,) 
x | (62) 
0 


tes 


¢C,, is the actual atom concentration at the wall and use has been made of the 

definition of a7,, from equation (41). p,,A,, can be found from equation (30). 
From the previous results the temperature variation along the wall can be 
found, using the thermal equation of state (equation (11)) for the mixture. The 
dimensionless temperature increment is given to a sufficient degree of accuracy by 
Cw — Cx 


63 
1+c,, 


As an illustration figures 4 and 5 show (c,,—c,,) and (7,,,—T.,,) (= 7',,7,,) plotted 
against x/7,, U fora @ of 5°. The values of p,, and 7., are taken to be one atmosphere 
and 4250 °K, and the gas is an ‘oxygen-like’ ideal dissociating one. In that case 
D’ in equation (30) is equal to 59,000/7., = 13-68 and with pz = 150g/ml., 
c, = 0-78. Thus p,,A,, = 0-086. 

The temperature rises steadily behind the corner (after the abrupt drop at the 
corner), approaching the asymptotic value shown in figure 5. The rise in tem- 
perature is a direct result of the fall in atom concentration. The recombination 
reaction results in the liberation of dissociation energy which reappears in the 
gas as random, thermal energy. 
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The energy equation (equation (14)) in the small disturbance approximation 
shows that A 4U°C,,,. (64) 


(N.B. C,,, = pressure coefficient at the wall.) Thus the enthalpy on the wall rises 
continuously after its abrupt fall at the corner and, in conformity with the 
adverse pressure gradient, the velocity must decrease as 2/7,,U increases (see 
equation (16)). The velocity on the wall behind the corner is always greater than U, 
however. 
Asymptote 
-0-02 


| —901 


0 04 08 1-2 16 20 
FicureE 4. Variation of atom mass fraction on the wall; B? = 1-5, (a,,,/a,..)? = 1-35, 
Po = latm., = 4250 °K, c,, = 0-780. = characteristics solution. 


— 200 
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Asymptote 
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Figure 5. Variation of temperature on the wall; B? = 1-5, (a,,,/a,,.)? = 1:35, 
Po = latm., T, = 4250 °K, c,, = 0-780. O = characteristics solution. 


From the values of p, 7’ and c computed above, the value of 7,, from equa- 
tion (39) can be calculated. In the present case it is found that 7,, = 3-85 x 109/k,... 
(It is found that 7’, of which 7,, is the mean value, varies by about + 10% for 
0 < 2a/7,,U < 2.) Since k, may be anywhere in the range 10'4-10" ml.?/mole?see, 7,, 
may be anywhere in the range 40sec to 0-04sec. For the example quoted 
U ~ 3x 10°cm/sec, so that 2/7,,U = 1 corresponds to an 2 of between 12 to 
0-012 em. The length scale of the relaxation zone is seen to depend quite critically 
therefore on the value of k,,.. Thus if the density increment p/,, could be measured 
(with an interferometer, for example) it should be possible to obtain a reasonable 
estimate of the magnitude of k,.,. 

Conditions like those adopted to illustrate the present theory could be achieved 
in the zone of equilibrium flow behind the primary shock wave in a shock tube. 
The feasibility of checking the theory by these means is being investigated. 

The linearizations which have been made in §3 would appear to be justified by 
the results obtained in the present section, provided both Mach numbers, M, and 
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M,, are neither too large nor too near unity. However, as 7,, becomes very small 
the gradients of velocity, temperature, etc., become very large in the streamwise 
(roughly, the x-wise) direction and it seems highly probable that it would be 
necessary in such cases to take account of the transport phenomena in the very 
short relaxation zone behind the corner. This question is, however, outside the 
scope of the present work. It should be reiterated that no account has been taken 
of the energy interchanges which occur in the internal modes of the molecules. 
The effects of these internal relaxations are ordinarily assumed to be negligible 
in the region of appreciable dissociation, but the present simple theory may be 
modified to take account of their presence and so provide a check on their relative 
importance. 


5. Further remarks on the flow round a corner 


The previous section has been almost entirely devoted to a discussion of the 
variations of p, p, 7’, etc., on the wall. This is primarily because it is much easier to 
evaluate these quantities analytically there, but solutions for the whole flow 
field can be obtained. In particular, the inversion of the exponential term in the 
transform of + and p can be accomplished (Morrison 1956). Reference to 
Morrison’s paper will show that the final (physical plane) values for the distur- 
bance variables would, however, be in the form of very intractable expressions. 
It is hoped that the simpler discussion of flow at (and hence, as can be seen from 
the form of the equations, also near) the wall will be of some interest. 

Finally, we may note some general results from the solutions obtained. The 
Laplace inversion theorem shows that 


1 (atin @ 
v'(,) E 2B,( dz. 


2m 


(ais areal number greater than the real part of the singularities in the integrand.) 
Closing the straight line contour with an infinite semicircle to the right, it can be 
seen that integration along the semicircle is equivalent to that on the original 
contour. The index of the exponential on the new contour tends to 


|2| exp [(targ z) 


as |z| > 00. Since — 47 < argz < 3m, the integral approaches zero, therefore, if 
£—B,n < 0. 

It follows that the effect of the corner is first felt along a line through the corner 
inclined at an angle sin~! (./;') to the free stream (i.e. the ‘frozen’ Mach angle). 
In contrast to the inert (or equilibrium) cases, it also follows from the integral 
above that the disturbance is not constant among characteristic lines passing 
through the corner (Kirkwood & Wood 1957, have shown that the characteristic 
directions are defined by the local values of .M, in a reacting system). This result 
is in line with the general result, quoted by the above-named authors, that simple 
wave flow does not exist in a reacting gas mixture. A possible physical explana- 
tion of this phenomenon in two-dimensional flow has been given after equation (56) 


above. 
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In order to investigate the flow in the vicinity of the equilibrium Mach line 
through the corner we may proceed as follows. Let us consider the pressure 
coefficient, which, from equations (52) and (53), is given by 


t l+o dw 


where we have written = w, = and = 9’. Then it is convenient to put 
= £'-7'B, =6 (66) 


and to rewrite the integral in equation (65) as 
20 ~ N \B?+0 P\B l+o 


x exp an —. (67) 


In this form a solution valid for large values of £’ and small values of 6 can be 
| obtained by the method of steepest descents. The index of the second exponential 
term in equation (67) is not a particularly manageable function of w but fortu- 
nately it turns out that w = 0 is a saddle point. The appropriate steepest path 
comes from Rl w = — © below the real axis, approaches and leaves w = 0 along an 
approximately parabolic curve given by Imw = + (—R1w)!and proceeds towards 
Rlw = —oabove the real axis. It includes a semicircular identation to the right 
of the simple pole at w = 0 and can be reconciled with the original inversion con- 
tour. Then we can write 
w{1 — = — 48%, 
where s is a real quantity, and it follows that 
(68) 
The term in braces in this integral must be expanded as a power series in the 


real variable s. A solution which is physically significant can be obtained from 
equation (68) when it is recognized that the expansion of the term 


+)/(1+o)} 
in the second exponential begins (isB?)/,/(B?—1). Then the leading term of the 
integral in equation (68) can be written as 


1 
exp (— }&’s®) sin {sdB/,/(B? — 1)}s-1ds 


1 
= — / 2 =~ 
The final solution for C,, can be written out in terms of 2, y, 7,,, etc., and is 


20 x—yB, ) 
3 + erf bez B-); (69) 
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Recalling the definitions of £’, T, ete., equation (69) may be expected to give an tl 

indication of the variation of C,, on either side of the equilibrium Mach line, si 

provided that x > 7,,U. (The writer is indebted to Prof. Lighthill for suggesting 

this approach to the problem.) re 
It should also be noted that the contour in equation (65) can be deformed into i) 

a single circuit embracing the branch points at w = —1 and — B?, plus a single z 

circuit surrounding the pole at w = 0. Then it follows that the latter circuit gives p 

rise to the dominant contribution (namely B-'), since any contribution from the p 

former must be O (exp (—&’)), when 9’ is appreciably smaller than £'/B,, and this st 

is negligible when &’ > 1. de 
It follows from equation (69) that at distances from the corner large compared re 

with the characteristic chemical length 7,,U, a smooth but quite rapid fall of 

pressure occurs across the equilibrium Mach line. The pressure coefficient would 

appear to reach its ultimate, equilibrium, value of — 24/6, without falling 

markedly below this, in contrast to the situation occurring on and near the wall. 

There, it will be recalled, the pressure drops discontinuously across the frozen Mach 


line to its frozen flow value, C,, = — 20/B,, thereafter rising steadily until the final 
equilibrium value is reached some way downstream. Equation (69) suggests that 
less and less pressure drop occurs across the frozen Mach line as distance from the 
corner increases, the majority of the decrease arising in the region centred about 
the equilibrium Mach line. The sharpness of this latter pressure drop for any 
given x and y is increased by a reduction of 7,, and it seems that the passage 
towards the full equilibrium flow limit of 7,,>0 occurs quite smoothly. The 


centring of the pressure drop about the equilibrium line will, however, occur for Al 
any value of 7,,, provided one is far enough away from the corner. | 
It is interesting to note that the dissipative effect of the chemical reactions in the te: 
gas smears the pressure drop over a wider and wider region around 2 = yB, as to 
x increases, under any given conditions. Thus, far from the corner, the flow is 
appears to be expanding through a fan of waves, even without the non-linear 7 
flattening effects which would be present in practice. One may conclude that A, 
such non-linear, convective terms will act to flatten out the pressure drop still ™ 
further in these regions. The present analysis is the two-dimensional analogue of - 
Chu’s (1957) treatment of the one-dimensional unsteady piston problem and the | W® 
conclusions are in agreement with those found in that case. ne 
mi 
6. Comparison with characteristics solutions 
Since the appearance of the first version of the present work (Clarke 19585) 
some characteristics solutions of the flow round a corner have become available | 
(Cleaver 1959). The results of this investigation are shownon figures 2 to 5forcom- | 
parison with the predictions of the linear theory. The characteristics solution was At 
A: carried out for an identical set of free-stream conditions to those used in the | 
linear theory example, and in fact constitutes an exact numerical solution of the 
equations 3 to 7 inclusive in § 2 above for the ‘oxygen-like’ ideal dissociating gas. | suf 
k, was assumed to be a constant, however, since no reliable information exists at act 
present of its temperature (and possibly also concentration) dependence. In be | 
view of the relatively small percentage changes in 7’ and c in the present example, par 


5 
} 
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this is perhaps not a serious objection, particularly since k, is unlikely to be a very 
strong function of either variable. 

The agreement between the exact and linear theories is seen to be quite 
reasonable, particularly when it is recalled that the linear theory will in any case 
overestimate the pressure drop at the corner even when the gas is chemically inert. 
The variation of pressure in the relaxation zone is predicted quite well by the 
present linear theory, as can be seen from figure 6 which shows the ratio p,,/P,,0 
plotted against 2/7, U for the two solutions (p,, = pressure on the wall, p,,. = pres- 
sure on the wall immediately behind the corner). With allowances for the classical 
defects in the linear theory, it would appear to give an adequate description of the 
relaxation zone behind the corner. 


Asymptote 


4 j i 
0 0-4 0-8 12 16 20 


FicuRE 6. Ratio of pressure on the wall p,, to pressure immediately behind the 
corner Py». O = characteristics solution. 


Appendix A 

The term 7~{A(1—c)—c?} which appears in equations (5), (6) and (7) in the 
text expresses the mass rate of production of atoms per unit mass of mixture due 
to the chemical reaction described in equation (1). It is derived as follows. The 
net rate of production of atoms measured in moles per unit volume per unit time 
is equal to 2{k, (concentration of A,) (concentration of A;)—k, (concentration of 
A,)? (concentration of A,)}, the concentration being measured here in moles per 
unit volume. To write this result in terms of mass fractions, it is observed that the 
mass fraction of A, = 1—c, and the mass fraction of A, = 1. The molecular 
weight of A, is the mean molecular weight of the mixture, namely W,/1+c. The 
net rate of atom production measured in terms of mass of atoms per unit mass of 
mixture per unit time is, therefore. 


W | p p 4p" 
— \k(l—c)— 1 
(4k W, 
W?2) \4k, p 
At equilibrium there is no net rate of atom production, and so 
(Wk,/4pk,), = c¢/1—¢,, 
suffix e implying equilibrium. Choosing ¢, to be the equilibrium value at the 
actual local p and 7’ values and making use of equation (11) to relate p, to p, it can 
be seen that k,W,/4pk, reduces to K as given above (equation (9)). The term in 


parentheses on the right side of equation (A 1) is 7~! as defined in equation (8). 
38 Fluid Mech, 7 
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594 J. F. Clarke 
The entropy equation (equation (7)) is derived from equations (14) and (15) 
The last term in parentheses in equation (A2) can be replaced by 
T1{K(1—c)—c?} 
by equation (5) and it can be shown that 
— fla = (R/W) log {c?(1 — c?)/c2(1 —c?)}; 
whence the form of equation (7) follows. 
The continuity equation in the present case can be written (see equation (57 a)) 


fou ov’ 
A3 
Plan ay) 
Writing density as a function of p, s and ¢, 
dp = dp + (5) 4 
cp os p,e 


s+ de. (A4) 
(0p/0p),.- = 47, the frozen sound speed. ne (A 


ce 


> 


) and (A3) we can now write 


4 2 On 

Oy Ox Cy (A5) 

Making use of equations (5) and (7), it follows that the last term in (A5) can be 


expressed as 


The quantity 7-'(0p/¢s), — T(es/ec),, ,} in (A6) is written as —po in 
equation (6) in the text. The value of o given in equation (10) can be obtained 
from manipulation of the thermodynamic equation. It is then in a suitable form 
for evaluation from the thermal and caloric equations of state. Details of the 
derivation are given by Clarke (1958a) but the equation of continuity was first 
given in this form by Kirkwood & Wood (1957). 


Appendix B 


The derivation of equation (58) from equation (57) proceeds as follows. Noting 


that 
we can write equation (57) as 
Pao pw lz pw (B 1) 
Thus Pa = 5 Cw + + Be ‘oxp(- (| dé". (B2) 
21 0 [ I 


Putting the value of C,,,, from equation (55) into the integral in equation (B2) 


pw 


gives a term shea to 


r *) exp - op 


+ | exp( | exp| | I, (B3) 
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Integrating the last term above by parts yields the terms 


2 
x I, | ae’ (B4) 
The last integral in expression (B4) cancels with the first integral in expression 
(B3). It follows from equation (B2) and equation (55) that 

2 2 


B, \ 


+B exp| - | | jaw. 
(B5) 


w 


The result quoted in equation (58) follows from reorganization of the terms in 
equation (B5). 

Equation (59a) for the concentration increment can be reduced in a similar way 
to that outlined above. Thus, from equation (55) we can write the last term in 
equation (59a) as 


9 B; 


exp [ —3(B?+1)u] Vuldua). 


0 0 


Integration of the last integral by parts shows that 


exp(—g/P) | aw 
~0 
20 


B,J, 


exp[—}(B?+ 1)W]/[3( 82-1) 
and the result quoted in equation (59) follows at once. 
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An experimental study of the forces generated by the 
collapse of transient cavities in water 


By I. R. JONES* AnD D. H. EDWARDS 
Department of Physics, University College, Aberystwyth 


(Received 20 July 1959) 


The paper describes an experimental investigation of the pressures developed 
at the seat of collapse of cavities in water. Single transient cavities, generated 
by a spark discharge, are allowed to collapse on the end of a piezoelectric 
pressure-bar gauge which measures the variation of thrust with time. It is 
shown that both the peak force and duration of the cavity collapse pulse are 
functions of the cavity lifetime. From an estimate of the minimum radius 
attained by the cavity and the peak force, the peak pressure on collapse is found 
to be at least 10,000 atm. Streak schlieren photographs of the collapse process 
show that a shock wave is radiated into the water at the moment of collapse 
and that the cavity rebounds. At the collapse of the rebound cavities the 
pressures developed are comparable with those developed by the collapse of the 
initial cavity, and these probably contribute materially to cavitational surface 
damage. 


1. Introduction 

Cavitation erosion was first observed by marine and hydraulic engineers at the 
turn of the century, and its study is of particular importance in these fields. 
During the last fifty years there has been an impressive accumulation of data 
and ideas concerning the mechanism of cavitation damage. A good deal of 
evidence, both theoretical and experimental, has accrued which supports the 
view that the origin of the surface damage lies in the large-amplitude brief 
pressure pulses which are generated when transient cavities collapse against the 
surface. The actual process by which metal is removed from the surface is not 
fully understood, but it may be conjectured to arise from an association of 
chemical, electrochemical, thermoelectric and mechanical effects brought about 
by these large pressure pulses accompanying cavity collapse. In view of the 
fundamental importance of the collapse pressure pulse in the damage process, 
there is a need for a detailed examination of its magnitude and character. 

A certain amount of experimental work has been reported on the observation 
of the pressure wave radiated in the liquid at the collapse of single transient 
cavities. Giith (1956) has obtained spark schlieren photographs of the radiated 
pressure wave, while Osborne (1947), Chesterman (1952), Harrison (1952) and 
Mellen (1956) have employed various types of piezoelectric pressure gauges to 


* Present address: Hydrodynamics Laboratory, California Institute of Technology, 
Pasadena. 


n 
2 
a 
: 
Vv 
| 
tl 
is 
al 
tl 
al 
at 
of 
th 
| | 2, 
| 
| 
| 
ele 
: | lo 
iy ca 
: du 
th 
an 
we 
we 
loy 
ba 
the 
| lea 
gai 
spi 
sin 
alo 
of 
du 
| 


Forces generated by the collapse of transient cavities 597 


measure the variation of the pressure, whose peak value is in the range 
2-10 atm, in the liquid at various distances from the seat of collapse of the 
cavity. It is doubtful whether such an experimental observation can furnish 
a reasonably accurate estimate of the peak pressure generated at the actual 
seat of collapse of the cavity. A more direct approach to the problem was pro- 
vided by the work of Ellis (1956) who demonstrated the existence in solids of 
strain waves originated by bubble collapse. Using an ultra-high-speed motion 
picture camera, ultrasonic cavitation bubbles were photographed collapsing on 
the surface of a photoelastic specimen, and photographs of the resulting transient 
isochromatic pattern in the specimen due to strains caused by cavitation were 
also obtained. A quantitative interpretation of the phenomenon indicated that 
the local stresses due to cavitation were of the order of 2 x 10° lb. in.-?, and a 
recording of the transient strains by a photocell showed that their duration was 
about 1 or 

The present study is concerned with the measurement of the forces generated 
at the collapse of single transient cavities on to a solid surface. The magnitudes 
of these forces were obtained by measuring the stresses developed in the solid 
during the coilapse of the cavity. In addition, streak schlieren photographs of 
the collapse phenomenon have been obtained using a rotating-drum camera. 


2. Experimental method 
2.1 The pressure bar 


A cylindrical pressure bar is employed to measure the forces developed by the 
collapsing cavities. Cavities are created on one end of the bar, termed the 
‘pressure-end’, by discharging a condenser, charged to a high voltage, through 
a gap between a tungsten needle and the end of the bar, which forms the earthed 
electrode. The energy liberated by the spark provides a source of extremely high 
local temperature in the region directly above the pressure end of the bar, thereby 
causing the water in this region to vaporize. Various sizes of cavity were pro- 
duced by varying the spark gap width, charging potential and the capacity of 
the condenser. Stress waves are propagated in the bar, both during the growth 
and collapse of the cavity, which can be detected and measured in a variety of 
ways (see Davies 1956). In the present investigation, resistance strain gauges 
were first employed, but were found to be unsuitable owing to their comparatively 
low sensitivity. This method was therefore abandoned in favour of the pressure- 
bar gauge described by Edwards (1958), which employs a quartz disk to measure 
the average stress over the cross-section of a composite bar of duralumin and 
lead. One difficulty which was encountered at first, both when using strain 
gauges and quartz disks, was the disturbing effect of the radiation from the 
spark on the gauge signal. This effect could not be removed by shielding alone, 
since the greater part of the disturbance was introduced by direct conduction 
along the metal pressure bar. In order to overcome this difficulty a short length 
of insulator, of the same diameter as the pressure bar, was inserted between the 
duralumin section and the quartz disk. To ensure undistorted transmission of 
the stress pulses through this insulator, its acoustic impedance must match that 
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of the duralumin. A heavy silicate flint glass was found to satisfy this require- 
ment; this glass was supplied by Messrs Chance-Pilkington (No. EDF 38919, 
Type 653/335). 

A diagram of the assembled pressure bar and its housing is shown in figure 1. 
‘Weston’ oil seals are used for mechanical support of the bar and also to provide 
a water-tight connexion between the bar and the tank. 
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FicurE 1. Diagram of the experimental tank, the spark-gap adjustment and the 
pressure bar and its housing. 


Pressure bars of two diameters, } and } in., have been used for measurement. 
These were calibrated by allowing shock waves in air of known strength, 
generated in a conventional shock tube, to be reflected at the pressure end of 
the bar. In addition, the dynamic behaviour of the bars was investigated by 
generating short duration stress pulses by the impact of steel balls on the 
pressure end. The observations are in agreement with those found by previous 
investigators (e.g. Ripperger 1952) that 4 and } in. diameter bars will faithfully 
transmit stress pulses of duration about 12 and 6ys, respectively. 


2.2. The tank and spark unit 


The experimental tank is constructed of 3°; in. thick sheet brass and is of square 
section of side 12 in. in length and 9} in. in height. Two circular water-tight 
observation windows, of 5 in. diameter, are set into opposite sides for purposes 
of photographic investigation; details are shown in the diagram of figure 1. 
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The open end of the tank can be covered with a sheet brass lid which is held 
securely in position by eight bolts. A device for varying the width of the spark 
gap, fixed centrally in the lid and insulated from it, enables vertical adjustment 
of the tungsten sparking needle, of diameter 0-08 cm, to be made. Vernier scales 
enable the width of the spark gap to be set to an accuracy of 0-01 mm. The 
spark is initiated by discharging a condenser by means of two thyratrons con- 
nected back-to-back (Mullard Type XH16-200); the design of the circuit is 
similar to the one used by Mellen (1956). Although oscillatory discharge is not 
prevented, the total discharge time of about 10s is sufficiently small, which is 
desirable in these experiments. 


2.3. The recording apparatus 


The voltage signal from the pressure-bar gauge is fed to the recording gear by 
means of a coaxial line contained in an earthed copper pipe. A wide-band 
amplifier, with a good high-frequency response, is used to amplify the signal 
before it is displayed on one beam of a double-beam cathode-ray oscilloscope; 
at the same time a timing wave is fed to the other beam from an oscillator. The 
time-base unit is initiated by a small voltage pulse which is derived from the 
spark-unit thyratrons when they become conducting. Since the interval between 
the voltage pulses obtained from the cavity growth and collapse is long com- 
pared with the duration of the pulses themselves, details of the recorded pulses 
would be lost if both were recorded on a long-duration sweep of the oscilloscope. 
Consequently, the time-base is designed to provide two successive sweeps when 
it is triggered, separated by a variable interval which can be determined 
accurately. In this way, the growth and collapse pulses are recorded on separate 
sweeps, while at the same time the time interval between them can be measured. 
A stationary plate camera is used to photograph the traces. 


2.4. The schlieren system and drum camera 


The schlieren system used for observing the shock wave radiated into the liquid 
when a cavity collapses is of the twin-lens parallel-beam type. The tank is placed 
between the two schlieren lenses, which have a focal length of 36 in. and aperture 
6in., with the glass observation windows aligned normal to the light beam. 
Variations of refractive index arising from density variations in the liquid are 
recorded as ‘streak’ pictures on a rotating drum camera. For this purpose a 
mask is placed over the tank windows in which a horizontal slit, 4 in. long and 
ig in. wide, is cut. The optical system gives an image reduction of }, and the 
maximum obtainable effective writing speed of the camera is about 0-7 mm/ys. 
The schlieren light source is a mercury flash lamp described by Edwards & 
Owen (1957) which is triggered by the spark-gap thyratrons. 


3. Results 3.1. The pressure measurements 
(a) The 4 in. diameter pressure bar 


In all the experiments described in this paper, the single cavities were created 
in tap water which had previously been allowed to stand in the tank for at least 
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24h to permit the gas content to reach equilibrium. The water used for the 
experiments was thus assumed to be saturated with air at room temperature 
and pressure. The results under these conditions were repeatable. 

The pressures recorded by the pressure bar are the mean values obtaining 
over the cross-section of the bar. Knowledge of the cavity diameter at any 


Spark gap Sparking Cavity 
width Capacity voltage lifetime Peak force 
(mm) (uF) (V) (us) (dynes x 10-*) 
0-01 0-01 1500 208 2-30 
0-02 0-01 2000 281 5-79 
0-02 0-25 1800 472 16-9 
0-01 0-50 2000 737 47-8 
0-02 0-50 2000 900 65-5 
0-04 1-00 2000 1059 838-6 
0-03 1:50 1800 1027 88-0 
0-04 4-0 1500 1463 105-9 


TaBLE 1. Typical results obtained for peak force and cavity lifetime for various 
sparking voltages, capacitors and spark-gap widths 
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Figure 2. Curves showing the relationship between the measured peak force F,,, of the 


collapse pulse and the cavity lifetime. ©, Data obtained with } in. diam. bar; @, data 
obtained with | in. diam. bar; +, data obtained from cavity rebound experiments. 


Ficure 3. Force-time curves for the collapse pulses (corresponding to regions AB and AB’ 


on curves of figure 2) obtained with } and ! in. diameter pressure bars. —, } in. diam. bar; 
cavity lifetime: 7804s, F,, = 57 x 108 dynes; ——-, | in. diam. bar; cavity lifetime: 360s. 


= 9:9 x dynes. 
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instant is required to relate the measured pressure to the actual pressure in the 
cavity at that instant; a discussion of the estimated actual cavity pressures will 
be given later. Thus if P denotes the applied pressure, assumed to be uniformly 
distributed over a small area a of the pressure end of the bar, then the total 
thrust F in the bar is Pa, and the corresponding voltage developed by the 
charge liberated by the quartz crystal is given by d,, Pa/C, where C is the total 
shunt capacity in farads and d,, is the piezoelectric modulus, whose value is 
2:2x10-"C/Kg force. The main observations that can be made from the 
pressure-bar measurements are of the character of the force-time relationship 
during the collapse stage of the cavity and of the variation of the peak force with 
the cavity lifetime, which is defined as the interval between the growth and 
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FicurE 4. Force-time curve of collapse pulse corresponding to the region BC of curve of 
figure 2, obtained with } in. diameter pressure bar. Cavity lifetime: 750ys, F’,, = 22-4 x 108 
dynes. 


collapse pulses. Table 1 gives a typical set of observations obtained from records 
taken using the } in. diameter pressure bar, and figure 2 shows the manner in 
which the peak force varies with cavity lifetime. The shape of the collapse pulse 
for a cavity of lifetime 780s is given in figure 3. 


(b) The 14 in. diameter pressure bar 
In order to establish whether the distortion of the recorded pulses, due to 
dispersion and attenuation of the stress pulses in the bar, was pronounced, all 
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the results obtained with the } in. bar were repeated using a } in. bar. The peak 
force vs. cavity lifetime relationship for this series of experiments is plotted for 
comparison on the same diagram, figure 2, as those obtained for the } in. bar, 
and the force-time curves for cavities of lifetime 360 and 750ys are given in 
figures 3 and 4, respectively. 


(c) Cavity rebound observations 

Following the collapse of a transient cavity, another cycle of growth and 
collapse may subsequently occur; this phenomenon of cavity oscillation, termed 
‘rebound’, has been observed by numerous investigators, e.g. Knapp & Hollander 
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Fiaure 5. Observed relationship between the rebound cavity lifetime T’ and the initial 
cavity lifetime 7’. 

Fiaure 6. Variation of peak force F’,, due to collapse of the rebound cavity with cavity 
lifetime T’. 


(1948), Harrison (1952), Chesterman (1952), and Mellen (1956). Evidence of 
rebound has been found for the cavities generated in the present work, and a 
brief investigation of the phenomenon was undertaken. The manner in which the 
rebound cavity lifetime 7’, defined as the time between the first and second 
collapse pulses, and the peak force F’,, due to the collapse of the rebound cavity, 
vary with cavity lifetime 7' is shown in figures 5 and 6. The values of (F/,, 7”) 
are also plotted on the curves of figure 2; they are seen to be in excellent agree- 
ment with the values of (F,,, 7’). 


3.2. The streak schlieren photographs 


A portion of a streak schlieren photograph obtained when a cavity of lifetime 
800s collapses on the plane surface of a } in. diameter pressure bar is shown in 
figure 7 (a) (plate 1); the main features of this photograph are labelled in the line 
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drawing of figure 7(b). A small length of the bar was allowed to project into the 
field of view, and this has given rise to the slightly dark band of constant width 
which runs down the length of the photograph. The completely dark band 
superimposed on the pressure bar shadow is due to the cavity, and its width 
corresponds to the diameter of the cavity at any instant. The estimated value 
of the minimum radius attained by the cavity is 0-049 em, which is approxi- 
mately the same as the radius of the tip of the tungsten sparking needle 
(0-04 em) which, in the present arrangement, inevitably appears in the field of 
view. For this reason a good estimate of the minimum cavity radius is not 
possible, but the value of 0-04 cm may safely be assumed as an upper limit; in 


Cavity Velocity 
lifetime of shock wave 
(#8) (m/s) 
350 1409 
445 1630 
800 1473 
815 1421 
843 1581 


TABLE 2. Measured values of average shock wave velocity for various cavity lifetimes. 
Velocity of sound in water at 13 °C = 1441 m/s 


all probability the true value is less than this. At the point of collapse a spherical 
shock wave is radiated into the water. The resolution obtainable on the records, 
however, does not permit of definite conclusions to be drawn concerning either 
the fine structure of the shock wave or the variation in its velocity in the close 
vicinity of the cavity wall. Average values of the velocity over a distance of a 
few centimetres obtained in a few experiments are given in table 2. It is seen 
that these values are sonic or just a little greater whereas near the cavity wall 
the records show that much higher values obtain. 


4. Discussion 4.1. The peak force-lifetime curves 


The peak force-lifetime curves, figure 2, obtained for the two pressure bars show 
that the peak force, for a given ambient pressure, is a continuous function of the 
cavity lifetime and hence of its maximum volume. Two distinct regions can be 
identified on these curves: (a) the portions 4B, AB’ which are concave upward, 
and (b) the flattened portions BC and B’C’ where the peak force increases com- 
paratively slowly with increasing cavity lifetime. The transition points B, B’ are 
fairly well defined and occur at roughly 1100 and 460s, respectively. Since the 
attenuation and dispersion effects of the stress pulses are different for the two 
pressure bars, the close agreement observed for the two curves in the region AB 
strongly supports the view that a true measure of the value of the peak force is 
provided by both bars. The discrepancies and change in character of the curves 
in the regions BC and B'C’ are attributed to the fact that the maximum cavity 
diameter exceeds the pressure bar diameter; consequently the collapse of the 
cavity will not be uniform and it will probably divide into a number of smaller 
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cavities. This fact has been confirmed by rotating-drum camera photographs. 
Moreover, it will be noted that the character of the collapse pulse in this region 
differs from that of the collapse pulse corresponding to the region A B’ ; whereas 
pulses in the latter region are smooth (figure 3), pronounced irregular oscillations 
are found on the initial portions of the collapse pulses for cavity lifetime in the 
regions BC, BC” (figure 4). 


4.2. The collapse pulse 


In the present work the duration of the collapse pulse is defined as the time 
interval, measured on the force-time curve, during which F rises above one- 
tenth its maximum value F,,. The two collapse pulses shown in figure 3, corre- 
sponding to cavity lifetimes of 780 and 360ys, are of 10 and 7ys duration, 
respectively. It was already verified that the frequency response of the pressure 
bar is adequate to deal with pulses of this duration; consequently these values 
could be taken to be a true measure of the pulse durations. This appears to be 
the first reliable experimental estimate of the collapse pulse durations. More- 
over, it is seen from figure 3 that the duration of the collapse pulse increases 
with the cavity lifetime. In §1 it was noted that Ellis & Sutton (1956, 1957) 
found collapse pulse durations of about 2s, which is a considerably lower value 
than those found here. The cavities studied by these authors were generated by 
a 20 ke’s sound field and had lifetimes of the order of 304s. Furthermore, they 
were created in fairly well degassed water (Ellis 1959). From this observation 
and the present results it may be inferred that purely vaporous cavities will 
give rise to sharper pressure pulses than gaseous ones. 

The shape of the stress pulse propagated in the pressure bar is similar to that 
observed by Mellen (1956) in the body of the liquid. A possible qualitative 
description of the phenomena, during various phases of the collapse process, 
which give rise to this type of pressure pulse is as follows. Initially the cavity is 
assumed to contain gas and water vapour, so that as the cavity starts to collapse 
the water vapour condenses at a rate governed by the wall motion. No appre- 
ciable effect on the solid will be observed during this stage. As the wall velocity 
increases the pressure of the gas inside the cavity begins to rise; this phase 
corresponds to the slow initial rise of pressure observed. When the wall velocity 
becomes appreciable, large pressures and temperatures are generated in the gas 
inside the cavity; this stage is observed in the solid as a sudden steepening of the 
front of the stress pulse. During this stage the wall velocity can reach super- 
sonic values depending on the initial gas content of the cavity. This is followed 
by a rapid deceleration of the cavity wall, as the minimum volume is approached ; 
at this stage a large compression pulse or shock wave is radiated into the sur- 
rounding liquid. 


4.3. The peak value of the collapse pressure 


Before the peak value of the pressure generated by the collapsing cavity can be 
estimated from the peak force measured by the pressure bar, a knowledge of the 
minimum diameter attained by the cavity is required. Furthermore, the fol- 
lowing three assumptions are made concerning the final collapse phase of the 


| 
| 
| ‘ 
d 
Vv 
| 
: 
t} 


Forces generated by the collapse of transient cavities 605 


cavity: (i) the force acts over a circular area at the end surface of the bar; 
(ii) the cavities remain in contact with the bar during collapse; (iii) the pressure 
distribution at the end of the bar is of the form shown schemiutically in the 
diagram of figure 8(a) in which the pressure is uniform inside the cavity and 
zero everywhere in the liquid outside the cavity wall. A sequence of spark 
shadow photographs of the growth and collapse of a cavity confirms the 
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FiGURE 8. Schematic representation of the pressure-distance distributions on a solid surface 
due to the collapse of a hemispherical cavity. The ideal distribution assumed in the present 
calculations is shown in (a) and that for a gas or vapour cavity in (b) where the maximum 
pressure occurs at the cavity wall; (c) refers to a Rayleigh cavity in which the maximum 
pressure is at some distance into the liquid. 


validity of assumptions (i) and (ii). Such a series of photographs is shown in 
figure 9 (plate 2) for a cavity whose total lifetime is approximately 800js and 
the pressure bar diameter is } in. (The time intervals between successive photo- 
graphs are not accurately known in this experiment but this information is not 
essential for present purposes.) These photographs clearly show that the cavities 
do in fact remain sensibly hemispherical during collapse and, what is of more 
importance, that they remain in contact with the bar at a stage approaching 
the final collapse volume. On this evidence the authors feel justified in concluding 
that the cavity remains in contact with the bar throughout most of its history, 
which implies that the total lifetimes of these cavities are too short for buoyancy 
forces to have an appreciable effect. 

It is difficult to assess the usefulness of the idealized pressure distribution on 
the end of the bar proposed in (iii) above. In the absence of any theoretical 
distribution, however, it appears to be the only reasonable assumption that can 
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be made at present. During the collapse of a gas or vapour cavity the pressure 
distribution on a surface is likely to be of the type shown schematically in 
figure 8(b), where the maximum pressure occurs at the cavity wall and there- 
after the pressure decays gradually in the region outside the cavity and not 
abruptly, as envisaged in the ideal case. On the other hand, in a Rayleigh cavity 
some arbitrary constant pressure obtains within the cavity throughout the 
collapse process and the maximum pressure occurs at some distance outside the 
cavity wall. The pressure distributions appropriate to this case are shown in the 
sketch of figure 8(c). In all probability the pressure distribution which actually 
occurs during collapse in the present experiments is a mixture of the type of 
distributions shown in figure 8(b) and (c). Consequently the estimated values 
of the peak pressures from the peak force data must be regarded as tentative. 
The method does, however, give an accurate force-time curve over most of the 
duration of the collapse pulse. 

Another factor which increases the element of uncertainty in the calculated 
peak pressures is the inconclusive nature of the information available about the 
minimum diameter of the cavity. A precise estimate of this quantity is difficult 
to make, and the various determinations made by various investigators differ by 
several magnitudes. For example, Harrison (1952) observed on his photographs 
that a fog of about 0-02 in. diameter remained after the disappearance of the 
cavity, and he therefore assumed that this represented the minimum diameter 
attained by the cavity. On the other hand, Sutton (1957) postulated that the 
peak collapse pressure acted on an area having a diameter of 0-001 in., this being 
the diameter of the small pits appearing on a plastic specimen subjected to 
cavitation. In the present study the experimental arrangement only permitted 
an upper limit of 0-032 in. to be defined for the minimum diameter, since the 
tip of the tungsten sparking needle interfered with the field of view. On this 
assumption the peak stress developed locally at the pressure end of the bar 
during the collapse of a cavity of lifetime 800s is 104 atm. This value, of course, 
represents the lower limit of the peak pressure and although it agrees with that 
obtained by Sutton, if this author’s estimate of the minimum volume were to be 
taken, then a predicted peak pressure of 107 atm would result. The minimum 
volume may be a function of cavity lifetime, and this is a point which requires 
investigation, consequently such a high value may be unrealistic. It is reasonable 
to assume, however, that from the results of the present work together with 
previous evidence, that the peak pressures developed are probably in the neigh- 
bourhood of 10° atm. Such large pressures would be expected to cause mechanical 
damage to a surface since they are well in excess of the yield strength of any 
material; e.g. the yield strength of the best steel is ~10' atm. Furthermore. 
from the few results which have been obtained for rebound, it is evident that 
the peak pressures developed at the collapse of rebound cavities can, and 
probably do, contribute materially to the disruption of a solid surface. 


4.4. The schlieren photographs 


The variation of cavity radius with time, derived from the schlieren photographs. 
is given in figure 10; the origin of time ¢ has been arbitrarily chosen such that 
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t =0 when the radius RF of the cavity is the maximum radius R, attained by the 
part of the cavity appearing in the field of view of the optical system. For the 
purposes of comparison, figure 10 also includes a computed curve showing the 
time history of collapse of a hemispherical cavity having the same maximum 
radius Ry and collapsing according to Rayleigh’s theory. Due to the spherical 
symmetry of Rayleigh’s result, the theory can be applied to the collapse of a 
hemispherical cavity having a zero or arbitrarily constant internal pressure on a 
semi-infinite body in an incompressible non-viscous liquid provided frictional 
losses at the liquid-solid interface are neglected. The experimental cavities are, 
of course, not of the Rayleigh type so that a comparison of the behaviour of the 


2 


Cavity radius R (mm) 


0 50 100 150 200 250 300 3 
Time (ys) 


0 400 450 


FicuRE 10. Variation of radius of cavity with time derived from the photograph of 
figure 7. Cavity lifetime 8004s; maximum radius of cavity 0-49 em. Experimental points, 
O—O. ———, Rayleigh hemispherical cavity. 


two is not strictly justified. Although the shapes of the experimental and the 
computed curves agree reasonably well, it is seen that the total collapse time of 
the experimental cavity is about 10°% shorter than that of the Rayleigh cavity. 
This result is surprising since one would expect frictional effects to prolong the 
lifetime of the experimental cavity compared with that of the equivalent 
Rayleigh cavity, and in fact Rattray (1951) has shown that the collapse time of 
a cavity situated near a wall is increased by about 20°% compared with that of 
a cavity in the body of the liquid. The discrepancy between the observed and 
experimental values must, therefore, be attributed to error in estimating the 
maximum radius of the cavity from the records. It will be noted that the 
smallest value of the radius measured is 0-049 cm; neither the quality of the 
record nor the fact that the radius of the tungsten sparking needle is itself 
0-040 cm permits a persual of the (R, ¢) analysis further than this point. 

In §4.1 it was concluded that a cavity which just totally covers the pressure 
end of a bar must have a lifetime corresponding to that of the transition points 
Band B’ on the (F,,, t) curves of figure 2. The lifetime of the cavity illustrated in 
figure 7 is 800s, and the maximum diameter attained is 0-98 cm. The assump- 
tion of a linear relationship between lifetime and maximum diameter implies 
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that a cavity having a maximum diameter of 1-27 cm (i.e. a cavity just covering 
the end of a } in. diameter bar) would have a lifetime of 1140s. This value 
agrees closely with the value of 1100s. appropriate to the transition point B’ 
of the (F,,, t) curve for a } in. diameter bar and thus confirms the correctness of 
the above conclusion. 


4.5. Gas content of the cavities 


In the preceding discussion no mention has been made of the influence of the 
gas content of a cavity on the dynamics of the collapse process. It is evident, 
however, that the gas content will have an important bearing on the pheno- 
menon. For example, Eisenberg (1953) has suggested that cavity rebound occurs 
if the initial gas content exceeds a certain lower limit, and Osborne (1947) finds 
that the peak pressure, measured at a distance from the cavity, decreases with 
increasing initial gas content. One important question which arises is whether 
the initial gas content decreases the damage potential of the cavity, i.e. does the 
presence of gas in the cavity attenuate the collapse pressure measured at the 
seat of collapse? In the present experiments rebound is observed, which implies 
that the initial gas content must be appreciable as would be expected, since gas 
is evolved by electrolysis and by release from solution. In spite of this fact, 
however, the pressures measured at the seat of collapse exceed the elastic limit 
of the material of the pressure bar. It may thus be concluded that damage can 
be produced by the collapse of cavities having an initial air content less than 
a certain upper limit as well as by the collapse of purely vaporous cavities. 


5. Conclusions 

It is shown that the pressure-bar method is suitable for measuring the peak 
force developed by the collapse of a single transient cavity in a liquid and for 
studying the shape of the collapse force-time pulse. The peak force in the collapse 
pulse is found to be a function of the duration of the pulse and of the cavity life- 
time. The minimum radius of 0-040 em attained by the cavity, estimated from 
the streak photographs, yields a value of peak pressure at the seat of collapse of 
104 atm. Consideration of the uncertainties in this estimate, however, together 
with values proposed by other investigators, indicate that the pressures are 
probably higher than this value and are more likely to be ~ 10° atm. Rebound 
of the cavity is observed both with the pressure and photographic recordings, 
and the peak pressure generated on collapse of the rebound cavity is of the same 
magnitude as that developed at the collapse of the initial cavity. Preliminary 
work on the observation of the shock wave radiated from the point of collapse 
show that its velocity at a short distance away from its point of origin is slightly 
greater than sonic. Although the recordings do not permit a precise estimate of 
the shock wave velocity near the cavity wall to be made, it is reasonably certain 
that here it is well in excess of sonic velocity. 


The authors wish to thank Dr T. G. Jones for his assistance in obtaining the 
schlieren photographs. 
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Figure 7 (plate 1). (a) Portion of a streak schlieren photograph of the collapse and 
rebound of a cavity of lifetime 800js; (b) labelled drawing of (a). 
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FIGURE 9 (plate 2). Series of spark photographs of the growth and collapse of a cavity 
generated by a spark discharge in water. Total cavity lifetime is roughly 800ys, and the 
pressure bar diameter is } in. 
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Forces generated by the collapse of transient cavities 
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Curvature of attached shock waves in steady 
axially symmetric flow 
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An electronic computer has been employed to calculate the ratio between the 
initial radii of curvature of the attached shock wave and the body for an axially 
symmetrical body in a uniform supersonic stream. The results are obtained with 
4 exact digits for more than 200 cases. They extend results obtained previously 
(Cabannes 1951) by means of numerical integration. 


1. Introduction 


We consider a body of revolution placed in a compressible fluid. The fluid 
possesses at infinity a uniform supersonic velocity ¢ parallel to the axis of revolu- 
tion Ox. A shock wave is formed in front of the body, and limits the region in 


Shock 


FicurE 1. Diagram of the flow field. 


which the flow is uniform. Viscosity and thermal conductivity are neglected 
outside the shock. We suppose that the surface of the body is tangential at the 
axis of revolution to a cone with semi-angle @,, and that the angle @, and the Mach 
number of the upstream flow have been chosen in such a way that the shock 
wave is attached at the vertex O of the obstacle. We locate the position of a 
point P in a meridian plane by the polar co-ordinates OP = r and 2 POx = 6 
(see figure 1). By means of these co-ordinates, the equation of the obstacle in the 
neighbourhood of the point O can be written in the form (1) and the equation of 
the shock wave, in the neighbourhood of the same point, in the form (2), namely, 
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The angle 6,, is determined by the theory of axially symmetric flow (Kopal 
1947); it depends on the Mach number M and the angle @,. The object of the 
present paper is to give tables for the determination of the value of the ratio 
(R/A) of the radii of curvature, at the axis of revolution, of the shock wave and 
the body; this ratio likewise depends on the Mach number M and the angle @,. 


2. Equations of motion 


We designate by wu and v the components of the fluid velocity at a point P in the 
directions 4 and (6 + 37), by p and p the pressure and density at this point, and by 
y the ratio of the specific heats of the fluid. The four functions w, v, p and p of the 
variables r and 6 satisfy the following partial differential equations which express 
the fundamental law of dynamics, the conservation of mass and the conservation 
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We attempt to satisty the preceding equations by means of functions expanded 
in series of whole and increasing powers of r, the coefficients depending only on the 
variable 0: ‘ 


v(r,0) = +5 + 
0 R | (4) 
p(r, 9) = +5 p,(0) 


r 


9) = pel) + 


By substitution of these expansions into equations (3) and by identification 
according to successive powers of r, one obtains an infinite set of differential 
equations. Equations (3) have a first integral, Bernoulli’s equation. As the 
limiting speed q,,, is constant in front of the shock and continuous across the shock 
wave, we have, valid in all the fluid, 

(5) 
y—Ip Im: e 


We also introduce the function a,(9) defined by: 
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The differential equations deduced from equations (3) can be written in the 
following form. Using given initial conditions, the functions with suffix 0 can be 
calculated from equations (7), while the functions with suffix 1 can be calculated 
from equations (8): 


Po (1 cot wht 
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3. Boundary conditions on the body 


The body is formed by the stream surface extended from the point 0. By 
expressing the condition that the differential equation of the stream function, 


(9) 


is satisfied by the function (1), one obtains the conditions 


v9(9,) = 0, (10a) 
v4(9,) Uo(9,) 
10 
* R ~ (106) 
According to the second of equations (7), one has that 7(9,) = — 2u9(@,): there- 
fore the condition (106) can be written in the form 
_ (11) 
R~ 3u(0,) 


4. Conditions on the shock wave 


At the shock wave, a certain number of conditions must be satisfied. These 
conditions, which express the fundamental law of dynamics, the conservation of 
mass and the conservation of energy, are expressed by equations (12), in which ¢, 
p and p designate the speed of sound, pressure and density in front of the shock 
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while £ is the angle which the tangent to the shock wave makes with the axis of 
revolution. .@ designates the Mach number along the normal (.@ = M sin f). 


~ 


—qsin@ a-4) cos — 0) 


(12) 
p yri y+ 
ytl 
The Mach number J is expressed as a function of the speed q by 
a2 
(13) 


By expressing the fact that the equations (12) are satisfied identically on the 


shock wave, one obtains the following values for the functions with suffix 0 and 1 


y—1¢@ cos?6 
= 
2 
p 2 I y-1 


tan = 0, 


2v,+ vy cot = 9, 
Y Ut 15 

Py - cot 0, ( >) 
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5. Numerical integration 

The numerical integration of equations (7) and (8) has been performed with the 
help of electronic computer gamma of the Faculty of Sciences of Grenoble. 'The 
great capacity of the machine and its high velocity of execution have allowed the 
computation of 209 cases to be performed, corresponding to 15 different bodies. 
The method of integration adopted is the Runge-Kutta method of fourth order, 
with intervals equal to one-twentieth of a degree; it seems that the value of the 
ratio of the curvatures can then be predicted with 4 exact digits. The results, 
which are given in the following tables,* have been computed with the adiabatic 


* For u,(9s)/Im = (1/6)! = 0-4082, the speed on the body, at the vertex, is sonic. 
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index having the value y = 1-4. The ratio of the curvatures is negative for the 
limiting velocity for which the shock wave is detached from the body; it is zero 
for a given value of the Mach number, which has been computed. 
In the case where the angle @, is small, it can be verified that the asymptotic 
formula, given by Rao (1956). 
R 40 1 
(16) 
MB 
is satisfactory for finite values of the Mach number. For higher values of @,, the 
results are exhibited graphically in figure 2. 


y= 14 
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2? 
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Ficure 2. Values of R/¥# vs M for various 6,. 
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Im M Im M R/2 
0, = 12:5° (cont.) 
tic 0-35 1-1732 89-224 —0-3234 0-4082 61-593 7-3049 
0-39 10215 86-921 —1-4832 0:45 1-2892 52-804 14-2637 
0-3913 11-0180 86441 —1-7700 0-5 1:4633 44-923 15-6587 
84-104 —2-2779 0:55 1-6623 38-934 13-1465 
16) 0-399 —-1-0168 80-841 6-0933 0-6 1-8916 34-144 10-0314 
0-4 1-0187 80-072 10-4070 0-65 2-1618 30-176 7-4305 
04082 1-0414 74-131 167-7417 07 2-4900 26-802 5:4927 
the 0-55 15151 41:363 4759-00 0-75 2-9070 23-869 4-0938 
0-6 17258 35-482 2196-71 0-8 3-4725 21-265 3-0785 
0-65 1-9699 30597 1172-21 0-85 4-3239 18-903 2-3215 
0-7 2-2611 26-372 576-76 0-9 5-8910 16-709 1:7300 
0°75 2-6224 22-592 255-98 0-95 11-1397 14-606 1-2335 
0:8 3-0961 19-107 103-76 
0-85 3-7723 15-792 29-586 0, = 15° 
0-9 4-8895 12-527 12-9795 0:3 1-2830 84-715 —1-0269 
0-95 7-4586 9-136 4-3398 0°35 1-1196 75-231  —0-8898 
0-99 22-6254 5-992 1:3343 0:3670 11289 70-070 (00000 
0. = 0-385 1-1608 64-965 1-6265 
0-4 1:1964 61-214 3-1447 
0-36 1-0980 87°324 0:6790 0-45 1:3448 51-391 6-6419 
0-6 17715 34-729 145-74 0-85 4-615] 20-587 1-7922 
0-7 2-3244 26-080 46-690 
0:95 16-8844 16-787 1-0793 
0:8 3-1985 19-326 13-516 0, = 17-5° 
0-85 3:9177 16-348 73241 03 1-2665 82-113 
0-9 51329 13-505 39809 0-3582 1:1707 68-912 —0-1426 
0-95 8/1036 10-665 2-104 03611 1:1743 68153  —-0-0000 
6, = 10° 0-45 1-4062 50-476 3-8682 
03765 1-0535 77°311 —0-8455 0-5 15881 44-067 4-2223 
0-3810 10576  75:347 0-0000 0-55 1:7996 39-089 3°9217 
0-4082 1-1190 64-899 15-9125 0-65 2-3438 31-756 2-9338 
0-45 1:2398 54713 42-7685 0-7 -28-939 2-4909 
0-5 14109 45-976 49-7619 0-75 3-1947 26-502 2:1100 
0-55 16049 39-445 39-7939 0-8 24-355 17717 
0-65 22-0872 29-915 18-1499 0-9 7-5793 20-699 12427 
0-7 2-4009 26-220 11-8199 6, = 20° 
75 9.7 9. 
ion. | 0-3 1-2728 79-220. — 11-1975 
eof | 0-9 54526 14-966 93434 0-3561 1-2259 66-670 0-0000 
0-95 91471 12-541 1-5017 
‘ 2 +345 56- *845 
0-45 11-4737 «50-008 2-6155 
= 12-5° 0-5 16609 44-205 2-8858 
0-3 1:3168 86-723 —0-8426 0:55 1-8806 39-658 2-7762 
0:3738  1-0902 72-409 0-0249 0-6 2-1404 35-980 2-5168 
0-4 1-1429 64-103 5-1224 0°65 2-4555 32-928 2-2497 
TABLE 1 


* The angles 4, are given in degrees. 
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M 
6, = 20° (cont.) 
2-8531 30-339 
3°3862 28-105 
4-1742 26-141 
5-5710 24-391 
9-6230 22-807 

= 22-6° 
1-3014 76-482 
1-2840 65-588 
1-2856 65-596 
1-3888 56-913 
1-4157 55-363 
1-5479 49-920 
1-7418 44-643 
1-9718 40-476 
2-2471 37-092 
2-5860 34-279 
3:0229 31-894 
3-6274 29-837 
4-5707 28-036 
16-7401 25-002 

25° 
1-3484 74-181 
1-3487 64-913 
1-3493 64-841 
1-3504 64-704 
1-3534 64-550 
1-4653 56-367 
1-4931 55-066 
1-6299 50-146 
1-8325 45-328 
2-0757 41-497 
2°3712 38-374 


2-7422 
+2335 
-9442 
-1439 
8-1165 
27-9654 
0, — 
1-4857 
1-5058 
1-5180 
1:6446 
1-6750 
1-8259 
2-0546 
2-3390 
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3°1745 
3-8645 
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24-0388 
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30-018 
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30° 
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63-068 
56°477 
55-409 
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37°202 
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34-171 
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69-681 
64-361 
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1-3591 
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— 1-1329 
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18415 
1-6881 
1-4817 
1-3309 
1-:0461 


— 1-:0138 
—0:0089 
0-0000 
0:0174 
0:0427 
1-0250 
11591 
1-5576 
1-7465 
1-7606 
1-6849 
1-5724 
1-4482 
13211 
1-2101 
1-0877 
0-8890 


— 00-7675 
00-0000 
0-1194 
0:7645 
0-8509 
1-1179 
1-2721 
1-3128 
1-2929 
1-2410 
1-1786 
1-1085 
1-0384 
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— 0:0047 
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0-3395 
0-35 
0-4 
0-4082 
0-45 
0-5 
0-55 
0-6 
0-65 
0-7 
0-75 
0-78 


0:3 
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0-34 
0-35 
0-4 
0-4082 
0-45 
0-5 
0-55 
0-6 
0-65 
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0-3383 
0-34 
0-35 
0-4 
0-4082 
0-45 
0-5 
0-55 
0-6 
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0-3 
0-32 
0-3395 
0-3396 
0-35 
0-4 
0-4082 
0-45 
0-5 
0-54 


0-3 
0-35 
0-4 
0-4082 
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6, = 35° (cont.) 
1-7105 64-316 
1-7309 63-009 
1-8766 57-548 
1-9115 56-656 
2-0875 53-216 
2-3633 49-758 
2-7228 46-957 
3°2112 44-653 
3-9340 42-729 
5:2059 41-122 
8-6893 39-706 
24-7546 38-964 
0, = 40° 
1-9533 69-453 
1-9938 65-177 
1-9982 65-147 
1-9993 64-902 
2-0235 63-863 
2-2027 59-282 
2-2463 58-527 
2-4712 55-606 
2-8441 52-647 
3°3752 50-238 
4-2147 48-299 
58293 46-589 
12-8630 45-186 
27-3166 44-907 
0, = 45° 
2-3720 70-058 
2-3910 68-077 
2-4336 66-381 
2-4387 66-228 
2-4718 65-354 
2-7233 61-487 
2:7866 60-847 
3°-1279 58-360 
3°7626 55-827 
4-8844 53-757 
7°7892 52-046 
26-2600 51-075 
50? 
3-1392 71-261 
3°1735 69-593 
3°2539 68-078 
3-2546 68-067 
3°3154 67-302 
3°7916 64-038 
3-9220 63-497 
4-7239 61-385 
70019 59-225 
21-4992 57-782 
= 65" 
5-4816 72-901 
6°2273 69-588 
9-8838 66-850 
12-0445 66°394 
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0-0000 
0-1325 
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0-6745 
0°8793 
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0-0133 
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0-8455 
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0:9217 
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The flow produced by an infinite rotating disk when the fluid at infinity is in a 
state of solid rotation is investigated numerically. When the fluid at infinity is 
rotating in the same sense as the disk, physically acceptable solutions exist in all 
cases. When the fluid at infinity.is rotating in the opposite sense to that of the disk, 
the only physically acceptable solutions appear to be those in which there is a 
uniform suction present acting through the disk. 


1. Introduction 

The steady flow of an incompressible viscous liquid, due to an infinite rotating 
disk, was first discussed by von Karman (1921). The liquid occupies the semi- 
infinite region on one side of the disk and the motion is rotationally symmetric. 
The effect of the disk is to throw the fluid near its surface radially outwards, and 
this in turn induces an axial inflow. The main interest in this problem is that, by 
virtue of assumptions about the velocity components, the Navier-Stokes equa- 
tions reduce to a set of ordinary, non-linear differential equations in a single 
independent variable. 

These equations are, in fact, the boundary-layer equations for the problem, since 
the terms which are ordinarily omitted in boundary-layer theory vanish identi- 
cally. Numerical integration of this set of equations thus yields an exact solution 
of the Navier-Stokes equations. Von Karman obtained an approximate solution 
to the problem using the integral method he invented, while Cochran (1934) 
corrected his solution and then calculated more accurate values by numerical 
integration of the equations. Bédewadt (1940) solved numerically the related 
problem of the flow produced over an infinite stationary plane in fluid which is 
rotating with uniform angular velocity at an infinite distance from the plane. 
Both these flows are particular cases of a general family of rotationally symmetric 
flows which were described qualitatively by Batchelor (1951). In these flows the 
fluid at infinity has an arbitrary uniform angular velocity about the axis of rota- 
tion of the disk, leading to a one-parameter family of solutions. 

Stewartson (1953) has also considered this problem and has obtained approxi- 
mate solutions in certain cases, while Squire (1953) linearized the equations for 


+ Present address: Atomic Energy Establishment, Winfrith, Dorset. 
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perturbations about a state of solid rotation. A complete discussion of these and 
other related papers is to be found in a review article by Moore (1956). 

Recently Fettis (1955) has devised an iterative process to obtain approximate 
solutions for several cases in which the fluid at infinity is rotating in the same 
sense as the disk. However, when the fluid at infinity is rotating in the opposite 
sense to that of the disk, the process diverges. It was therefore decided to investi- 
gate numerically the general set of equations, and to find out for what cases 
physically acceptable solutions exist. All the numerical integrations have been 
performed on the Pegasus digital computer at the University of Southampton. 

It should perhaps be emphasized that the present work, like the above- 
mentioned papers, refers to an infinite disk and is further based on the assumption 
of a similarity solution. Consequently, any edge effects for a disk of finite radius 
are automatically ignored, and in cases when the flow in the boundary layer is 
radially inwards the question of the range of validity of the similarity solution has 
still to be settled; a discussion of this point has been given by Stewartson (1957). 
When the flow in the boundary layer is radially outwards the similarity solution 
is simply the first approximation in an expansion of the dependent variables in 
powers of r. 

Section 2 introduces the notation and derives the equations. §§3 and 4 contain 
certain asymptotic and approximate solutions. The numerical methods employed 
are discussed in $5, and in §6 the necessary extensions are given so that suction 
through the disk may be included. The results are given in the concluding sections. 


2. The equations governing the motion 


Cylindrical polar co-ordinates (r, ¢, z) are used, with the disk in the plane z = 0, 
and the fluid occupies the region z > 0. Assuming the similarity solution, with the 
dependent variables in the form 


u=rQF(é), v=rQG(0), w= (vQ)H(6), p/p = vOP(C)+ 4xQ?r?, (1) 


where x is a constant and ¢ = z(Q/v)}, we find (see, for example, Schlichting 1955, 
p. 75) that the Navier-Stokes equations become 


F" = F?-G?+HF' +k, (2) 
= 2FG + HC’, (3) 
= (4) 
2F +H’ =0, (5) 


where a dash denotes differentiation with respect to €. Equations (2), (3) and (5) 
must first be solved for the three velocity components F, G and H, and then the 
pressure distribution can be found immediately since (4) integrates to give 


P, = —}3H?. 


The boundary conditions which must be satisfied at the disk are 


u=0, v=rQ, w=0 
and, in view of the similarity assumption (1), these are equivalent to 
F(0)=0, G(0)=1, A(0)=0. (6) 
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At infinity, the fluid is rotating with uniform angular velocity sQ and hence 
F(x) =0, G(o)=s. (7) 


However, these conditions are not sufficient to determine the solution uniquely, 
for the constant « in equation (2) is still arbitrary. It is also necessary to assume 
that both F’ and F” vanish at infinity, and thus it follows from (2) that 


K == (8) 


It will be seen, in fact, that the asymptotic solution which satisfies (7) auto- 
matically makes the velocity derivatives vanish at infinity, thus ruling out any 
possibility of a shear layer in the fluid at a great distance from the disk. Thus, the 
simplified set of Navier-Stokes equations to be solved is 


F" = F2?-@?+HF' (9) 
G” = 2FG+ HC’, (10) 
2F +H’ = 0, (11) 


subject to the boundary conditions (6) and (7). 

When s > 1 it was found more convenient to use the angular velocity of the 
fluid at infinity as a reference velocity. Denoting this by w, and the angular 
velocity of the disk by ow, a dimensionless co-ordinate is defined by £ = z(w/v)}. 
The velocity components are 


u=rwF(é), w= 


and the differential equations become 


F" = F?-GFt HF +1. (12) 

=2FGLHYGY, (13) 

2F + H' =0, (14) 

where dashes now denote differentiation with respect to £. The boundary condi- 

F(0)=0, G(0)=0, #(0) =0,) 
and =1. 


Since ow is the angular velocity of the disk, the range 0 < o < | corresponds to 
values of s greater than 1. In particular, the value ao = O gives the problem already 
solved by Bédewadt, and his solution shows that boundary-layer effects extend 
out to about € = 8. 

Various characteristics of the flow patterns can be computed, in particular the 
torque on the disk. Assuming for the present a finite radius R, and assuming the 
similarity solution is accurate over an appreciable part of the disk, we find that the 
moment is given by R 

M=- 


II 


—7R4(vQ)? G'(0). 


This quantity will be required subsequently in a discussion of the numerical 
solutions. 
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3. Asymptotic solutions 
For large values of z, the flow is very nearly uniform solid rotation, and an 
asymptotic solution can easily be found. If we write 


= f(¢), G = s+g(C), H = h(¢), (16) 


where c is the component of the velocity at 0 in the axial direction towards the 
disk, and is a function of s, and then neglect squares and products of f, g and h and 
their derivatives, equations (9), (10) and (11) become 


f" = —289—<f", (17) 
g” = +2sf—cg’, (18) 
= 0. (19) 
It follows from these equations that 
f = cos Bsin nC} (20) 
and g = {Boos A sin (21) 
4 + + 
where A= 3|-e- + | | (22) 
1 ((c4 + — c2)4 


It should be noted that A is always negative, irrespective of the sign of ¢ and the 
magnitude of s. A and B are arbitrary constants whose values will be determined 
later for various values of s by fitting the numerical solution. 

These expressions display the oscillatory nature of the flow away from the 
disk and (23) shows the dependence of the wavelength on both the axial flow and 
the angular velocity of the fluid. These solutions are the leading terms in formal 
expansions for F and G in powers of e’¢. These expansions have been calculated by 
Bédewadt up to and including terms of order e®¢, This was necessary in his work 
because the asymptotic solution was matched with a power series solution at 
€ = 1. However, in the present work the asymptotic solution is fitted to the 
machine solution at € = 6, and then only as a final check on the solution. Conse- 
quently, only the first term is required. Asymptotic solutions of (12) and (13) will 
also be needed later: the derivation of these is precisely similar to the above 
results, and it is found that 


F = ek{ cos Asin py, (24) 
and G = sin 25) 
l + 64)3 + 
where A, = 5 9 (26) 
1 + 64) 
and 5 (27) 


As before, —c is the limiting value of the axial component of velocity at infinity. 
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4. Approximate solutions 


It will be seen in the next section that, for the numerical method employed, 
estimates of the velocity derivatives at the disk are needed. One way of obtaining 
these, of course, is to use the von Karman integral method, but in practice it was 
found more accurate to linearize the equations about a state of solid rotation. 
When s or o take values close to unity, the fluid at infinity is rotating with nearly 
the same angular velocity as the disk. Writing o = 1 +6 there are formal expan- 
sions for the functions .F, GY and # as follows 


F =F +8F (28) 
G = (29) 
HK = KH 0H 4+ PH ot... (30) 


When 6 = 0, the state of solid rotation is recovered and thus 


F,=0, G=1, #,=0. (31) 
In bey of the boundary conditions for .*, Y and #, the boundary conditions on 
F,,.G, and #,, are 
F,(0) =F,(0) =0 (n=1,2,3,...), 
G (0)=G,(0)=0 (n= 2,3,...), 
= = ( ) 
G(0)=1, G,(w) =0, | 
H,(0) = 0 (n = 1,2,3,...). 
If the above expansions are now substituted into (12) to (14) and coefficients of 


powers of d are equated, sets of linear differential equations are obtained for the 
functions F,, G,, #,,. The first few solutions are found to be 


F, =e sin€, 


G, = eScosé, (33) 
H , = 
ee 
2 cos £ + sin £) cos 5, (34) 
25 
2 = (cos& + sin &) 3sin &), 
and for the third set 
+3 
Ft+tG,= ( (l+2 — (+ 
400 100 


(35) 
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The first approximation (33) was obtained by Squire (1953). It will be noticed 
that for large values of £, away from zeros of sin € and cos €, the dominant terms 
in F, and Y, are those involving €; similarly, in (35) the dominant term is of 
order £7. It can be shown that the dominant term in 7, +79, is 


— 1)! exp[—(1+72)&] 


and thus for large values of £, except near = nz and (n+4)7, where n takes 
positive integral values 
F+iG =~ i+dexp[—(14+ 46)é] {sin€ +7 cos (36) 


The importance of these expansions lies in the fact that they provide fairly 
accurate values of F'(0) and Y’(0) which are required in the numerical integration. 
From (33), (34) and (35) 


F'(0) = 6+ 0-16? + 0-012568 + O(44), (37) 
G'(0) = —d—0-26? + 0-022563 + O( 54), (38) 
= —8 + 0-362 — 0-0875683 + O(54). (39) 


Despite the fact that these expansions can only be expected to give reasonable 
accuracy for values of |6| which are small compared with unity, on putting é = — 1, 
so that the disk is at rest, these give the very accurate estimates 


F'(0) = —0-9125, G'(0) = 0-7775, = 1-3875, 
which may be compared with the final values, obtained numerically, namely 
F'(0) = —0-9420, G'(0) = 0:7729, = 13696. 


On the other hand, the form of equation (36) shows that trouble may be expected 
when 6 < —2 and this is borne out subsequently. 

A similar linearization can be carried out for equations (9) to (11); writing 
s = 1+¢e, we find that 


a = F'(0) = —e—0-4e? + 0-0625e3 + O(e4), (40) 
b = G'(0) = € + 0-3e? — 0-0475e63 + O(e4), (41) 
—c = H(x) = e—0-2¢? + 0-0125e3 + O(e?). (42) 


These are not so accurate as expressions (37), (38) and (39) but nevertheless 
provide suitable starting values for the subsequent calculations. 


5. Numerical method of solution 


The differential equations requiring numerical methods for their solution are 
those given in §2, (9)-(11). Solutions of these equations are required which 
satisfy the boundary conditions (6) and (7). For the purposes of computation on 
Pegasus, it is necessary to ensure that all quantities are less than unity in absolute 
value, so the following substitutions are made: 


02-4, = F2*, yy=H2-%, ys = P24, ye = G'2~. 
(43) 
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With these substitutions, the equations (9)-(11) may be rewritten as a system 
of six first-order differential equations. This reduction is performed because the 
standard Pegasus subroutine may then be used to solve them. The equations in 
question are: 


W=24 Ys= —Yo } 44 
Ys = — 277-8 + yyy, 2748729, = + 27. 
The first of these is used to keep track of how far the integration has proceeded. 


The associated boundary conditions are obtained by substitution of (43) into (6) 


and (7). Thus, 
0, y(0)=90, y,(0)= y,(0) = 0; 


= 0, = 82-4. 


These six conditions are sufficient to define the problem mathematically but it 
must be remembered that, in §§2 and 3, the additional restrictions that y;(0o) 
and y,(00) must be zero were imposed in order to obtain physically acceptable 
solutions. 

In the numerical work, the values of the scaling factors were chosen to be 
a=0, =2, y =d=3 and e = 4. These choices, which were suggested by the 
solutions given by Cochran (1934) and Bodewadt (1940), proved to be satisfactory 
for all values of the parameter s. 

Since the functions vary fairly rapidly for small values of € and more slowly for 
larger values, it was found convenient to use the Runge-Kutta method of 
integration as modified by Merson (1958). The modification suggested by Merson 
provides an automatic change of interval (either a decrease or an increase) when 
necessary. 

In order to start the integration of the system (44) at ¢ = 0, it is necessary to 
estimate the values of y;(0) and y,(0); these are denoted by a and b (see (40) and 
(41)). Reasonably good estimates can be made using these equations in the cases 
when s = +1. With these values, the system was integrated from zero to € = 12. 
This was repeated using a + Aa, 6, and a, b+ Ab, as assumed values for y;(0) and 
y,(0). (Aa and Ab are small changes in a and 6.) In this way, three sets of values of 
y,(12) and y3(12) were obtained. These were used to evaluate 


Ca” ch Ca cb 


Then, using the first-order terms in a Taylor expansion, it is possible to solve a pair 
of linear equations and to calculate da and db, which are corrections to be applied 
to the original a and 6 in order to make y,(12) and y3(12) more nearly 0 and s2-%, 
respectively. 

Such a sequence of operations is called an iteration in the sequel. The iteration 
was repeated until a and 6 were obtained to an accuracy of about ten decimal 
places; in most cases this made y,(12) zero to an accuracy of about + 10~* and 
likewise y,(12) differed from s2-4 by about + 10-*. 

This technique worked well fors = +0-9 ands = + 0-8, but unfortunately the 
series (40) and (41) did not produce sufficiently accurate initial values of a and 6 
for s < +0-8. (The solutions diverged before the integration reached ¢ = 12.) 
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The obvious way to obtain more accurate values is to take more terms in the 
series but the algebra involved in obtaining further terms theoretically was 
prohibitive. Consequently, the coefficients of e* and e° in (40) and (41) were found 
numerically using the final, accurate values of a and b for the cases s = +0-9 
(=e = —0-1) ands = +0-8(=6¢ = —0-2). The more accurate series obtained in 
this way are 

F'(0) = —0-4e2 + 0-0625¢5 — 0-0179764 + 0-0074865 + O(E°), ) 

G'(0) = +¢+ — 0-0475e3 + 0-0169964 — 0-00894e° + O(e8). J (4) 

Such a method enabled solutions to be obtained fors = +0-7 and + 0-6, but for 
s < +0-6 the initial values for a and b calculated from (45) allowed the exponen- 
tially increasing part of the solution to enter before the integration had reached 
€ = 12. Further coefficients in (45) could not be obtained using the previous 
method because all the significant figures cancelled and the resulting coefficients 
were very inaccurate. 

Thus, a different approach was required for s < + 0-6. Von Karmén’s solution 
was first improved by iterating with s = 0. Then the values of a and b so obtained 
were used to obtain an approximate solution in the case s = +0-1. As was to be 
expected, divergence occurred for quite small values of € (of the order of 4). So the 
boundary conditions y,(%) = 0 and y,;(%) = were satisfied at € = 3 by 
iterating as before. The revised values of a and 6 so obtained enabled the integra- 
tion to be extended to € = 4 before divergence set in. In this way, iteration was 
applied at € = 4 to make y,(4) = 0 and y,(4) = s2-?. This process was repeated 
until € = 12 was reached. At this stage, the changes in a and b from € = 11 to 
€ = 12 were O(10-*) and the calculation was stopped. This method was applied 
successfully for s = + 0-1 to +0-5so covering the range 0 < s < +1-0. 

Theranges > +1 was treated by considering the equations (12)—(15) written in 
terms of ¢ = 1/s. Obviously the original programme needed only trivial modifica- 
tion to do this. To obtain starting values for this set of equations, the series 
expansions (37) and (38) in terms of 6 (= o — 1) proved accurate enough to obtain 
solutions throughout the range 0 < 7 < +1 but two more terms of the series were 
calculated numerically as before. The results are: 


F'(0O) = +8+0-16? + — 0-0126364 — 0-0044665 + 
G'(0)} = | — § — 0-24? + 0-022568 + 0-00380d4 — 0-001156° + O(6*). (46) 
So far in this section, attention has been confined to positive values of s and a. 


In order to obtain solutions for —1 < s < 0 the same technique as that used for 
small positive s was employed. The values of a and b for s = 0 were used as initial 


estimates for s = — 0-05. This appeared to work satisfactorily for s = — 0-05, 
—0-10ands = —0-15, but fors = — 0-2 the values of F’(12) and G’(12) were very 
different from zero. Thus, this solution (for s = —0-2) was unacceptable on 


physical grounds.+. It was apparent that something peculiar was happening even 


+ The authors have been informed that in an unpublished investigation A. C. Browning 
found that the solution with s = —0-2 is unacceptable on physical grounds. Browning 
also considered other values of s or o and his solution for o = 0 is reproduced in 
Schlichting (1955). Our results for this case agree better with those of Bédewadt (1940) 
than with those of Browning. 
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before the solution was printed out, because in going from the iteration at € = 11 
to that at € = 12 the changes in a and b were not small. 

Despite the fact that this solution was unsatisfactory on physical grounds, the 
method was continued to s = — 1-0 at which stage the values of a and b were very 
small, O(10~4), but F’(12) and G’(12) were O(1). Thus the region of large 
gradients, or boundary layer, seems to have moved away from the disk and large 
shears are appearing within the body of the fluid at the end of the integration 
range, where there is no solid boundary. 

The solutions obtained for s = — 0-05, — 0-1 and — 0-15 appear satisfactory at 
first sight (that is, #(12) and G(12) are O(10~4)), but on closer examination they 
display several peculiar features. In the first place, it is shown in §2 that the 
torque on the disk is proportional to |@’(0)| = |b]. When s = 0, |b| = 0-615922 
and when s > 0 the value of |}| is less than 0-615922, which is to be expected on 
physical grounds. Now it is also to be expected that |6| will increase as s becomes 
negative, since the disk is being rotated in the opposite sense to that of the fluid at 
infinity. It appears from the numerical solution that |b| in fact decreases as 
s becomes negative; for example, when s = — 0-05, |b| = 0-615676, and when 
s= —01, |b| = 0-608253. The velocity components for these cases also attain 
their free-stream values at considerably greater distances from the disk than in 
the cases when s is positive or zero. It may be that the situation would be clarified 
if the range of integration was extended; on the other hand, the changes in a and b 
in going from the iteration at € = 11 to that at € = 12 were small. An alternative 
boundary condition at infinity is to postulate that the velocity derivatives vanish 
there, but in this case the constant « in equation (2) is no longer known, and 
problems concerning the uniqueness of the solution then arise. The flow between 
two rotating disks is being currently investigated and it is hoped that this will 
throw further light on the problem. 

It should perhaps be emphasized that great care was taken to ensure that a 
physically sensible solution had not been missed. For example, the value of b 
which was guessed in the case s = —0-05, was originally taken to be more 
negative than — 0-615676, but the ‘solution’ converged to the original one in all 
cases. 


6. The problem with suction 


Stuart (1954) has solved a similar problem with s = 0 and suction applied at 
the disk. As was to be expected, suction stabilizes the flow; that is, it reduces the 
magnitude of the radial and angular velocity components. The question thus 
naturally arises as to whether the application of suction to the disk when s is 
negative, might prevent the boundary layer moving off the disk. Suction through 
the disk can be applied, in the present notation, by making y,(0) negative instead 
of zero. Following the notation of Stuart, a suction parameter a* can be intro- 


duced such that 
H(0) = —a*. (47) 


The equations governing the motion are, as before, (9)—(11), but the boundary 
conditions are (6) modified by (47) and (7). 
40 Fluid Mech. 7 
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The value s = — 1 was chosen as a reasonably typical one and various values 
of the suction parameter a* were tried. It was found that the values 1-5, 1 and 
0-8 all gave physically acceptable solutions and quite clearly larger values of a*, 
corresponding to a stronger suction, would also be satisfactory. However, for 
a* = 0-4 no acceptable solution could be found: as before, the iterative process 
led to a solution for which the boundary layer is no longer attached to the disk, 
but appeared at the other end of the range of integration. In order to obtain a 
qualitative picture of the situation, the von Karman integral method was 
employed in the following way. Equations (9) and (10) were integrated from 
€ = 0 to € = 4, subject to the boundary conditions 


F(0)=0, G(0)=1, 


F(3)=0, G(d)=s, (48) 
F’(é) = 0, @'(é) = 0. 
On making use of equation (11), it follows that 
- F’(0) = | (3F2 —G?) d+ (49) 
~0 
G'(0) = 4 | | Fd€+a*. (50) 
J0 J0 
A quartic polynomial was assumed for F and a cubic for G, and on puttings = — 1, 
equations (49) and (50) reduce to 
—a 19a? 39 
= + 51 
~ 210 70’ 
3 13 
Ps (52) 


where a = dF’(0). 
These equations give imaginary values of « if a* < 0-257, but real solutions for 
all values of the suction parameter greater than this. For very large values of a*. 
(51) and (52) show that 


F’'(0) ~ 


G'(0) ~ — 2a*. (53) 


and these are in qualitative agreement with the values obtained by expanding 
F and @ in inverse powers of the suction parameter, as was done by Stuart for 
the case s = 0; the derivatives are 


F'(0) = 2 +0/ ): = (54) 


a* a* 


7. Results 


Many results have already been mentioned in $$5 and 6, but some graphs and 
tables are presented here to emphasize the more important points. Figures 1—4 
show the radial and transverse velocity profiles for different values of s and o. 
The velocity derivatives at the disk are given in Tables 1 and 2 for the indicated 
values of s and o, respectively. The constants describing the flow at a great 
distance from the disk are presented in Tables 3 and 4. These include the axial 
velocity of the fluid, the damping factors A and A, defined in equations (22) and 
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(26), the wave-numbers yw and , defined in (23) and (27), and the constants A, 
B, ¥ and &. It will be noticed that the values for A, B and ¢ for s = 0 differ 
slightly from those obtained by Cochran. 

An example of a flow in which suction at the disk is occurring and the fluid at 
infinity is rotating in the opposite direction to the disk has been computed; 
Ae 


oF 


0:10 0:20 F(¢) 
Radial velocity profiles 0 < s < 1 


Figure 1. The function F defined in equation (1), 0 < s < 1. The value of s is indicated 
in the figure. The radial velocity of the fluid is rOF. 


8 F’(0) =a =b 

0-0 +0-510233 — 0-615922 
0-1 + 0513397 — 0601554 
0-2 4.0:501870 — 0-572080 
0-4 + 0-439580 — 0-478673 
0-6 +0-331470 — 0348434 
0-8 40-183469 — 0187591 
0-9 + 0-095936 — 0-096951 
1-0 +0-0 —0-0 


TABLE 1. Values of the velocity derivatives F’(0) and G’(0) for various values of s. 
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complete numerical solutions for this case, and for the cases s = 0-6 and e 
o = 0-8 without suction, are being held at the Editor’s office for consultation by in 
interested people. 0 
8. Conclusions 
When the uniform rotation of the fluid at infinity is in the same sense as that of 
the disk, a physically acceptable solution exists for each value of s. In every case 
there is a boundary layer attached to the disk and in all cases, except one, the 
motion approaches the uniform rotation at infinity in an oscillatory fashion. The 
6 
Transverse velocity profiles 0 < s < 1 

FicurE 2. The function @ defined in equation (1); the zonal velocity of the fluid is rOG. Fre 
These are all cases when the angular velocity of the disk is greater than that of the fluid at 
infinity. 

F'\0) G'(0) 

0-0 — 0-941971 +0-772886 

0-1 — 0844923 + 0-718393 

0-2 — 0-751682 + 0:658418 

0-4 0-569080 + 0-522477 

0-6 — 0:385194 + 0:366431 

0-8 —0-196121 +0-191812 

0-9 — 0-099014 + 0-097977 

1-0 —0-0 +0-0 1 

TABLE 2. Values of the velocity derivatives F’(0) and Y’(0) for various values of o. 


at 


Flow in the presence of a rotating disk 629 


exception to this is the flow (for which s = 0) originally discussed by von Karman 
in which the fluid at infinity is not rotating but moving in a purely axial direction. 
One unexpected feature of the solution is that the axial velocity at infinity for 
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Radial velocity profiles 0 < 7 < 1 


Ficure 3. The function F defined in equation (12); the radial velocity of the fluid is 
rwF , and the appropriate values of o are indicated on the figure. 


8 c= —H(o) —A —p A B 

0-0 0-88446 0-88446 0-0 +0-91772 +1-20211 
0-1 0-91769 0-95931 0-19981 + 0-48981 + 1-32529 
0-2 0-86175 0-99062 0-35730 +0-19045 + 1-20306 
0-4 0-65996 1-00683 0-59097 + 0-01320 + 0-75996 
0-6 0-43077 1-00510 0-75977 — 0-00993 + 0-45192 
0-8 0-20801 1-00146 0-89141 — 0-00332 +0-21003 
0-9 0-10201 1-00037 0-94800 — 0-00261 +0-10431 
1-0 0-0 1-0 1-0 —90-0 +0-0 


TABLE 3. The flow at a great distance from the disk: these are the constants in (20) 
and (21) which give the approximate velocity distribution for large values of ¢. 
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8 = 0-1 is stronger than the corresponding value for s = 0. The maximum value of 
the radial velocity is also greater in this case. This feature is also displayed by the 
approximate results of Fettis. In the range 0-2 < s < 1, the flow relative to the 
disk gets progressively weaker as s increases to unity. No such anomalous 
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Transverse velocity profiles 0 < 7 < 1 


Ficure 4. The function Y defined in equation (12). 


—c = 4H(o) —A, — fly -B 
0-0 1-36961 0-43841 0-89031 0-23543 1-02911 
0-1 1-19516 0-49529 0-91502 0-18303 0-88037 
0-2 1-03080 0-55306 0-93593 0-13694 0-75559 
0-4 0-72608 0-67043 0-96761 0-06787 0-55281 
0-6 0-45378 0-78606 0-98721 0:02698 0-37316 
0-8 0:21272 0-89647 0-99718 0-00749 0-19283 
0-9 0-10309 0-94912 0-99934 0-00128 0-09882 
1-0 0-0 1-0 1-0 0-0 0-0 


TABLE 4. The flow at a great distance from the disk: these are the constants in (24) 
and (25) which give the approximate velocity distribution for large values of &. 
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behaviour is found in the range 1 > o > Oand this may have some bearing on the 
fact that the expansions in powers of 6 (equations (28)—(30)) are rather more 
accurate than the corresponding expansions in powers of ¢. In any case, the 
linearization technique is shown to give an accurate picture of the flow for quite 
a wide range of values of the expansion parameters 6 and e. 

When the fluid at infinity is rotating in the opposite sense to that of the disk, 
and s < —0-2, then it appears that no steady solution is possible unless there is 
suction acting, which prevents the boundary layer from leaving the disk and 
attaching itself to the ‘disturbing agency’ at infinity. For each value of s in this 
range there is probably a critical value of the suction parameter a*, giving the 
minimum amount of suction needed to keep the boundary layer on the disk. All 
that can be said is that for the case s = — 1, this critical value lies in the range 
0-8 > > O-4. 

The situation as regards values of s in the range 0 > s > — 0-2 is not clear: 
solutions of the equations satisfying the postulated boundary conditions appear 
to exist, but they give anomalous values for the torque on the disk. 


The authors would like to express their gratitude to Prof. Howarth for several 


helpful discussions. 
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A laminar planetary jet 


By ROBERT R. LONG 


School of Engineering, The Johns Hopkins University 
(Received 15 July 1959) 


A theory of a jet in a rotating, viscous fluid is developed. It is suggested that it 
may be related to a jet observed in an experiment with a rotating spherical shell 
of liquid and, in addition, may partly explain the existence of the subsurface 
equatorial current in the Pacific Ocean. The fundamental physical idea is that 
of a balance of vorticity brought into the jet by advection and diffused by 
friction. The necessary approximations are borrowed from boundary-layer 
theory. The linear case is solved completely and inferences about the non-linear 
jet are obtained by dimensional reasoning. 


1. Introduction 


One of the intriguing features of large-scale fiow patterns in the oceans is the 
tendency for some of the surface and subsurface currents to concentrate in narrow 
strips or jets. The best known of these is the Gulf Stream, but a few observations 
have been made recently of a westerly current of planetary scale in the Pacific 
several hundred feet below the surface at the Equator. It is called the equatorial 
undercurrent or Cromwell current and is described in a number of papers (Crom- 
well, Montgomery & Stroup (1954), Knauss & King (1958), Hidaka & Nagata 
(1958)). Briefly, it is an eastward moving current of water 300km wide sym- 
metrically located with respect to the equator. The maximum speed is about 
3 knots. Its longitudinal extent is still in doubt but it extends at least 3000 miles 
west of the Galapagos Islands. 

Velocity concentrations occur also in the atmosphere and in certain laboratory 
experiments. One of interest to this paper was found in an experiment by the 
author (Long, 1952) with a shell of water between two rigidly connected con- 
centric hemispheres in rotation. An obstacle is moved along a latitude circle in 
the water in the direction of the rotation but either slower or faster. If the 
obstacle rotates slower than the fluid, a series of waves develops around the globe 
as in figure 1. If the obstacle moves faster than the fluid, it drags around with it 
the westerly jet? of figure 2. 

A plausible explanation for the jet in figure 2 (and a parallel explanation for the 
absence of the jet in figure 1) may be developed from the schematic drawing in 
figure 3. We suppose that the motion in the shell is two-dimensional, i.e. parallel 
to spherical surfaces and with no variations in directions normal to these surfaces. 
Ahead of the obstacle the fluid moves more slowly than the barrier and speeds up 


+ This is really wake flow rather than jet flow, but the distinction is not important for 
our present purposes. 
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to a maximum just in the lee. Outside, fluid particles move slowly toward the jet 
from pole and equator. If we neglect lateral friction, except in the jet, particles 
bring in the absolute vorticity? of the latitude at which they originate. This means 
an influx of positive relative vorticity on the poleward side and negative relative 
vorticity on the equatorial side. This influx can support the strong shear at the 
edges against the diffusing effect of viscosity. 


N 


S 
FiGuRE 3. Schematic jet. 


2. Theory 


We approach this problem analytically by assuming a viscous fluid moving 
two-dimensionally in relative motion parallel to a spherical surface with a basic 
rotation 2. Density variations are zero. The radius of the sphere is a, latitude 0, 
longitude ¢, the potential of external forces V, density p, pressure p, constant 
viscosity v, and velocities 

dé dé 


u* =acosd—, v* =—. 


dt dt 


The Navier-Stokes equations of motion in the rotating co-ordinate system are 


acosO a acos@ acos0 
* * 
(1) 
azcos@ a®cos?@ a®cos?@ a?cos?@ %)’ 
u*vt sind +V 
5 +. 2.0u* sin = _(Ple+V)o 
acos@ a acos@ a 
(2) 
a a?cos@ a2cos?@ azcos?@ | 


+ We refer to the vorticity component perpendicular to the sphere. This is conserved in 
frictionless, two-dimensional flow. 
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The equation of continuity is 


ve 
— + = 
acosd acosé 
so that = —* u* = alee 
a COs a 


where 7* is the stream function. 

In accordance with the discussion of §1, we assume a jet flow moving west- 
east at latitude 4, with a width of the order 6. It is convenient to adopt two new 
independent variables «*, y*, where x* is distance along the axis of the jet in the 
direction of motion and y* is distance toward the north from latitude 0). Let 


dx* =acos0,d¢, dy* = add, 
and expand sin # and cos@ in a Taylor series about 4, 


y* cosd,y** 
cos 4 = cos4,—sin of 
a # 
y* sind, y* 
= sin + cos 4,” 
a 2 a 
In these expansions y*/a is very small} in the jet and its vicinity, and from the 
viewpoint of boundary-layer theory we should at first glance feel justified in 
writing cos? = cos4,,sin@ = sin@, everywhere in the equations of the problem. 
We will do this except in the Coriolis force term where we retain the first two terms 
in the expansion of sin#. This must be done because the first approximations to 


the Coriolis terms, — 2Qv* sin 4), 2Qu* sin 4, can be written 
act 


and are consequently eliminated with the pressure terms when we cross- 
differentiate to produce the vorticity equation. The vorticity effect of the earth’s 
rotation,{ discussed in $1, would therefore be lost unless we retained an addi- 
tional term in the expansion of sin @. 

Using this approximation and dropping other terms that are small from the 
viewpoint of boundary-layer theory,§ the equations become 


+ — = (3) 
ye j 

py*u* = — Xp. (4) 

+ 0, (5) 


+ This is certainly true in atmospheric and oceanic jets. The jet of figure 2 is relatively 
wide and extends around the globe. Our theory which ultimately reduces the globe to a 
plane, except for the dynamic effect of the variation of angular velocity, may have only 
a loose connexion with the experimental jet. 

t Rossby (1939) was the first to use this approximation for the Coriolis parameter 
20 sin 7. It has been used frequently in meteorology and oceanography, usually without 
attempts to justify it. 

§ It is an interesting, and perhaps an important fact, that the largest neglected terms in 
our approximate theory are one order smaller for an equatorial jet than for a jet located at 
an arbitrary latitude. 
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where = 20 cos @,)/aand = p/p + V —2Qsin 6) The kinematic momentum 
transfer is L 


If we require that w* and its first derivative with respect to y* approach zero as 
ly*| + co, J must be constant. We can now eliminate all constants from our 
problem by defining the non-dimensional quantities 


u* 
Ke ps 
—, £= 7 
where K = —J.The resulting equations have the same form as (3)—-(6) except that 
the constants are all missing: 
UU, + VU, — YU = —X,+Uy,, (8) 
yu = (9) 
+0, = 0, (10) 
—l= 2| (v—yyr + u*) dy. (11) 
0 


Thus the problem of the jet can be solved once and for all, and for all conditions, 
by solving equations (8)—(11). This is not done in this paper although we indicate 
below a sequence of successive approximations and solve completely the first 
approximation. At this point we draw some inferences about the Cromwell jet 
from order of magnitude arguments applied to (8)—(11). Since the boundary 
conditions imposed on the problem introduce no new constants, we may be 
confident that all quantities w, v, x, x, y, entering (8)—(11) are of order of magnitude 
one. Assuming this and applying equation (7) we evaluate certain quantities 
of interest in the special case of the equatorial undercurrent. If we take 
1-5 x 10?emsec™! as an observed speed at some point on the axis of the jet and 
B = 2Q/a ~ 2-3 x 10-3 em— sec, the first relation in (7) yields 


K ~ 5-7 x 10% see. 
Then 6 ~ 2:6 x 107em 


is an estimate of the width of the current. This is of the same order as the observed 
width and the numerical value should be compared with a calculation below using 
the linearized equations. It isremarkable that the above estimate does not require 
any knowledge of the assumed coefficient of friction. It can be verified from the 
linearized theory that the phenomenon as a whole is very insensitive to the 
magnitude of the friction.t A further check on the order of magnitude of 
quantities in the undercurrent is in the computation of the horizontal length 


scale from (7). We get 1017 
Lw- > em? sec}. 


ft One might also guess that the phenomenon is insensitive to the exact form of the law 
of friction. This would be important if true because the eddy viscosity concept, especially 
with a constant coefficient, is admittedly crude. 
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If we use a figure of 7 x 107 cm*sec™! for v, suggested by Montgomery & Palmén 
(1940), we find Z ~ 1-4 x 10®cm or perhaps 8000 miles. This seems at first glance 
a little large, by a factor of two or three, but not if we append the observation that 
measurements show little or no variation of the current over a distance of 
3000 miles or so. 

Equations (8)~(11) may be put in a better form by using as dependent variables 
y and A, the latter being defined by 


A, = (12) 
A(a, 0) = 0. 


Then the problem is determined by 


= yy t 


= (13) 
—A,,— = 0, (14) 
A(x, 00) = 1, (15) 


together with the conditions of symmetry and finiteness. 

We may approach the solution of (13)-(15) by noticing that far enough up- 
stream the velocities in the jet are arbitrarily small and only linear terms are 
important. The linear problem has solutions of the form 


Ao h,(z), 
= (—2x)-41,(2), 
where z=(-a)ty. 
We are led then to a scheme of successive approximations 
A= (2). (16) 
> (—2)-#2+5)] (2), (17) 
n=U 


Substituting (16) and (17) into (13) and (14), we obtain sets of equations in the 
single independent variable z. The first set is 


—hiz = 0. (18) 
ho — 21, = 0, 


where the prime denotes differentiation with respect to z. 


3. First approximation 

The system (18) may be solved very simply by expanding in a Taylor series 
about z = 0. After satisfying the symmetry and finiteness conditions at z = 0 
two constants remain. One may be determined in terms of the other by imposing 
the finiteness condition at z = 00. The single remaining constant can be found by 
satisfying h)(0o) = 1. In the numerical computation, one obtains sufficient 
accuracy by assuming /) = 0 at z = 10 or so, since imposition of the condition at 
both z = 10 and z = 8 does not give significantly different results except near the 
extreme value of z. The solution is shown in figure 4. 
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The u-component reveals counter currents on both sides of the main jet with 
maximum speed about one-half that of the primary jet. The motion dies off 
rapidly further north and south. The width of the primary jet is given approxi- 


mately by ~ (19) 
It narrows very gradually downstream. The width increases with the viscosity 


and decreases with higher rotation, but very insensitively. In the ocean at 
x* = 5x 108em, the width is 8 x 107 em, more than twice that of the undercurrent. 
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Ficure 4. Velocity profile in the linear portion of jet. Since we have not imposed the 
condition that hy (00) = 1, the abscissa scale may be multiplied by an arbitrary constant. 


4. Non-linearity 
The solution of (8)-(11) in regions where the non-linear terms are important 


has not yet been found. In the case of the equatorial undercurrent the non-linear 
terms will be unimportant if, for example, v*uj./fy*v* < 1, or 


Using u* ~ 102? emsec"!, ~ 10-8 ~ 2-6 x 10%em, we 
see that B ~ 1 and the condition is not met. The fact that the non-linear terms 
are not overriding (as they are, for example, in atmospheric jet phenomena) 
indicates that the ocean current may be qualitatively similar to that in the linear 
case. The effect of non-linearity can be judged by considering the physical picture 
of the maintenance of the jets by a balance of advection and diffusion of vorticity. 
We see that in the linear solution the jet system narrows as / increases. The non- 
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linear terms serve to bring in additional perturbed vorticity and their inclusion 
should have effects similar to those associated with an increase of £, namely, a 
narrowing of the jet. The modification is in the desired direction because the 
linear current turned out to be too wide to fit the observations of the undercurrent. 


The author wishes to thank Prof. R. B. Montgomery for several stimulating 
conversations about the equatorial current. The research was sponsored by the 
U.S. Weather Bureau and the Office of Naval Research. 
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